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Abstract 


Several  past  research  efforts  have  developed  methods  to  take  advantage  of  Global 
Positioning  System  (GPS)  positioning  and  apply  it  to  a  precision  landing  system  (PLS). 
There  have  been  proposals  to  phase  out  the  current  Instrument  Landing  System  (ILS)  in 
favor  of  a  more  cost-efficient  and  effective  system.  Accomplishments  have  been  made  in 
detailing  a  system  implementing  an  INS  aided  with  differential  GPS,  a  GPS  pseudolite, 
and  a  radar  altimeter  to  handle  the  critical  PLS  requirements.  This  research  applies  the 
newly  developed  Modified  Multiple  Model  Adaptive  Estimation  (M^AE)  architecture  in 
an  attempt  to  enhance  the  performance  of  a  PLS  in  an  environment  involving  GPS 
interference.  The  M^AE  uses  Multiple  Model  Adaptive  Estimation  (MMAE)  and 
Kalman  filtering  techniques  to  estimate  the  levels  of  interference  and  the  navigation 
performance  of  the  aircraft  simultaneously.  In  addition,  in  the  original  development  of 
M^AE,  the  truth  and  filter  model  used  were  of  the  same  order.  This  research  serves  as  a 
demonstration  of  M^AE  applied  to  system  where  the  truth  model  is  of  a  higher  order  than 
the  filter  model. 
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Integrated  GPS/INS  Precision  Approach  Landing 
with  M^AE  Interference  Avoidance 


1.  Introduction 


The  focus  of  this  research  is  to  examine  an  application  of  a  newly  developed  method 
of  simultaneous  state  and  parameter  estimation  to  the  precision  approach  and  landing 
problem  [1 1,  22, 45].  The  study  involves  a  simulated  aircraft  navigation  system  using  an 
inertial  navigation  system  (INS)  augmented  with  data  from  altimeters  and  Global 
Position  System  (GPS)  receivers.  It  extends  the  efforts  of  previous  analyses  from  the  Air 
Force  Institute  of  Technology  (AFTT)  entailing  integrated  aircraft  navigation  systems  [7, 
13, 34,  38, 56]. 

1.1.  Background 

There  has  been  extensive  research  in  aircraft  automatic  precision  approach  and 
landing  systems  [5, 7, 10, 11, 13, 18,  22, 45, 56].  These  systems  prove  invaluable  during 
periods  of  low  visibility,  or  when  a  pilot  is  unfamiliar  with  the  landing  site.  Landing 
safely  involves  many  factors,  including  flying  at  the  proper  speed  and  attitude,  knowing 
where  the  aircraft  and  the  landing  point  are  accurately,  and  having  a  clear  view  or 
indication  of  the  landing  point. 
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1.1.1.  Instrument  Landing  System 

The  current  system  in  use  by  the  Department  of  Defense  and  civil  aviation  is  the 
Instrument  Landing  System  (ILS)  [10].  This  system  is  a  highly  directional  guidance 
system  that  gives  the  pilot  critical  information  pertaining  to  the  aircraft’s  azimuth  relative 
to  the  centerline  of  the  runway,  called  the  localizer,  and  the  aircraft’s  approach  elevation 
relative  to  the  surface  of  the  runway,  called  the  glideslope.  The  modulated  signals 
provide  aircraft  course  and  direction  information  to  the  pilot  [47].  The  system  transmits 
two  principal  carrier  signals,  108. 10  MHz  to  111.95  MHz  for  the  localizer  and  329. 15 
MHz  to  335.00  MHz  for  the  glideslope.  The  localizer  beam  can  be  received  35°  from 
both  sides  of  the  runway  centerline  up  to  10  nautical  miles  (nmi)  away  and  10°  up  to  18 
nmi.  The  glideslope  beam  is  1.4°  wide  centered  at  an  ideal  glideslope  of  approximately 
3°.  This  translates  to  approximately  a  1,000  ft.  drop  in  altitude  for  every  3  nmi  traveled, 
a  guideline  many  pilots  use  for  a  quick  check  [58].  In  addition,  this  signal  is  normally 
usable  up  to  10  nmi  away.  Both  carrier  signals  are  modulated  with  a  different  frequency 
signal,  depending  on  which  side  of  the  beam  centerline  the  aircraft  is.  The  two 
modulating  frequencies  are  90  Hz  and  150  Hz.  (Figures  1  and  2  illustrate  how  this  is 

applied). 

This  system  is  in  use  extensively  and  satisfies  the  Federal  Aviation  Administration 
(FAA)  Category  IH  requirements  for  aircraft  precision  approaches  [12].  However,  due  to 
the  cost  of  maintaining  the  aging  ILS  structure,  there  is  an  active  search  for  replacement 
systems.  One  alternative  being  actively  pursued  involves  applying  Global  Positioning 
System  (GPS)  based  systems  as  a  low-cost,  flexible  substitute  [5, 18, 45]. 
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Figure  2. 


Glideslope  Service  Range 


3 


GPS  is  a  well-established,  reliable  system  for  navigation  and  surveying,  and  is  currently 
being  used  as  a  supplemental  navigation  aid  for  non-precision  aircraft  approaches. 

1.1.2.  Global  Positioning  System 

The  FAA  is  working  closely  with  state  aviation  and  industry  officials  to  develop 
GPS-type  procedures  for  all  U.S.  airports  [5, 18, 45].  Several  procedures,  unique  to  GPS, 
have  already  been  developed.  One  method  uses  differential  GPS  (DGPS)  to  provide  error 
correction  data  to  incoming  aircraft  that  use  GPS  as  a  navigational  aid  [5,  18, 45]. 

Another  advantage  of  GPS  is  its  ability  to  provide  the  pilot  more  flexible  landing 
approach  options,  unlike  the  ILS,  which  obliges  landing  aircraft  to  follow  one  prescribed 
approach  pattern  closely,  limiting  the  pilot’s  choices. 

Any  alternative  navigation  system  must  still  conform  to  the  FAA  precision  landing 
requirements  and  must  satisfy  Category  HI  performance  before  it  can  be  considered  as  a 
replacement  to  ILS.  Table  1  gives  the  constraint  parameters  for  the  different  categories 
[7, 13]. 


Table  1.  FAA  ILS  Precision  Approach  Requirements  (la) 


Category 

Horizontal  Accuracy 

Vertical  Accuracy 

I 

±  28.1  ft 

±  6.8  ft 

II 

±  8.6  ft 

±  2.8  ft 

III 

±6.8  ft 

±1.0  ft 
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Relying  solely  on  a  GPS-based  system  is  inadequate.  The  GPS  transmissions,  which 
are  universally  available,  are  those  of  the  standard  positioning  service  (SPS).  As  a 
security  measure,  SPS  can  be  degraded  by  introducing  a  random  drift  known  as  selective 
availability  (SA)  to  reduce  the  precision  available  to  users  outside  of  the  U.  S.  Armed 
Forces  [37].  Currently  SA  is  active,  and,  as  shown  in  Table  2,  the  level  of  accuracy  for 
GPS  is  unacceptable  for  even  a  Category  I  approach.  However,  even  without  SA,  the 
precision  available  fails  to  meet  any  of  the  FAA  criteria  for  vertical  accuracy  (see  Table 
2).  This  lack  of  precision,  along  with  the  unavailability  of  the  precise  positioning  service 
(PPS)  signal  to  the  public  has  encouraged  the  development  of  other  methods  to  decrease 
the  error  levels  in  SPS.  One  method  that  had  gained  popularity  is  the  use  of  differential 
GPS  (DGPS). 

DGPS  yields  greater  position  accuracy  by  eliminating  errors  common  to  GPS 
receivers.  The  reference  DGPS  receiver  is  located  at  a  surveyed  spot,  known  with  a  high 
degree  of  accuracy,  e.g.,  at  an  airport.  When  the  GPS  signals  it  receives  are. compared  to 
those  from  another  receiver,  e.g.,  on  board  an  aircraft,  error  corrections  can  be  calculated 
and  sent  to  the  other  receiver,  improving  the  accuracy  of  its  positioning;  see  Figure  3. 
The  two  receivers  should  be  within  27  nmi  (50  Km)  to  be  effective  [24].  While  highly 
accurate,  positioning  using  DGPS  falls  short  of  the  required  precision  for  FAA 
approaches.  Table  2. 

Centimeter-level  precision  can  be  achieved  using  another  method  known  as  earner- 
phase  GPS,  which  calculates  the  phase  shift  of  the  SPS  carrier  signal,  LI,  at  1,575.42 
MHz.  The  wavelength  of  LI  is  approximately  19  cm;  the  positioning  accuracy  is  a 
fraction  of  this  wavelength.  However,  determining  the  number  of  whole  (integer)  cycles 
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between  the  satellite  and  receiver  becomes  a  problem.  A  momentary  loss  of  the  signal 
can  disrupt  the  count  of  the  whole  cycles,  resulting  in  integer  ambiguities  [4].  Currently, 
most  carrier-phase  uses  have  concentrated  on  surveying  applications,  since  determining 
integer  ambiguities  is  more  difficult  in  a  dynamic  environment  [41, 42]. 


Table  2.  Precision  Values  Available  with  GPS 


GPS  Method 

Horizontal  Accuracy 

Vertical  Accuracy 

SPS  with  SA 

±  134.8  ft 

±  168.6  ft 

SPS  without  SA 

±  33.5  ft 

±42.0  ft 

SPS  using  DGPS 

±  7.2  ft 

±  9.2  ft 

Reference 

Receiver 


Figure  3. 


Differential  GPS 


1.1.3.  Inertial  Navigation  System 

Most  aircraft  use  an  inertial  navigational  system  (INS)  for  their  standard  navigation 
aid.  An  INS  is  a  self-contained  “box”  (it  requires  no  outside  information)  that  uses  a 
system  of  gyros  and  accelerometers  to  sense  specific  forces  as  the  aircraft’s  position, 
velocity,  and  attitude  changes.  With  a  knowledge  of  the  gravitational  field,  the  INS 
integrates  the  specific  forces  over  time  to  calculate  its  velocity  and  position  with  respect 
to  a  known  reference  frame.  A  significant  drawback  of  an  INS  is  that  the  readings  in  the 
vertical  direction  (relative  to  the  earth)  are  unstable.  This  is  usually  corrected  by  the 
addition  of  a  barometric  altimeter  to  the  INS.  For  better  accuracy  at  lower  altitudes,  a 
radar  altimeter  can  further  improve  the  vertical  position  accuracy.  Another  problem  is 
that  INS-derived  positions  have  errors  that  grow  over  time,  resulting  in  a  drift  [27].  This 
can  be  corrected  using  ground-based  navigation  aids  such  as  VHF  Omnidirectional 
Range  (VOR)  and  Tactical  Air  Navigation  (TACAN). 

GPS  can  serve  as  a  replacement  for  these  ground-based  aids  and  adds  flexibility  in  the 
flight  path,  as  the  aircraft  does  not  have  to  follow  a  trajectory  to  stay  in  range  of  a  VOR 
or  TACAN  system.  As  seen  above,  GPS  also  has  a  weakness  in  vertical  positioning. 

This  is  due  to  the  fact  that  the  local  horizon  of  a  receiver  limits  the  visibility  of  the  GPS 
constellation.  Only  satellites  above  the  horizon,  or  more  commonly,  above  a  user- 
selectable  elevation  mask  angle  (set  to  minimize  ground  clutter  and  atmospheric 
conditions)  can  be  used  to  determine  a  navigation  solution.  One  method  to  improve 
vertical  GPS  precision  uses  ground-based  GPS  transmitters,  called  pseudolites,  to  provide 
a  GPS-like  signal  from  below  an  aircraft’s  horizon.  Pseudolites  require  extremely  low 
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power  transmissions  to  mimic  the  —163  dBw  GPS  signal  power  level  while  not  saturating 
the  LI  frequency  with  its  own  signal  and  blocking  actual  GPS  satellite  transmissions. 

The  low  power  level  of  GPS  transmissions  brings  to  light  the  problems  of 
interference  and  noise.  GPS  spread  spectrum  provides  some  protection,  still  some 
interference  has  been  experienced  from  users  when  receiving  spurious  emissions  near  the 
LI  band,  such  as  second  and  third  harmonics  of  UHF  television  stations,  and  even  higher 
order  harmonics  from  VOR  and  TACAN  [40].  Also,  there  is  the  risk  of  someone 
deliberately  broadcasting  interference  to  deny  any  GPS  signals  to  a  receiver  (“jamming”) 
or  mimicking  transmissions  from  a  satellite  to  give  erroneous  information  (“spoofing”). 

Recall  that  an  INS  is  a  self-contained  system  that  does  not  rely  on  outside 
navigational  signals  or  beacon  transmissions.  Despite  its  inherent  long-term  drift 
characteristics,  an  INS  provides  velocity  and  position  estimates  with  high  degree  of 
precision  and  low  error  standard  deviation  in  the  short-term.  A  great  deal  of  work  has 
been  done  integrating  systems  such  as  GPS  receivers  and  one  or  more  INSs  (for  failure 
redundancy),  combining  the  benefits  of  GPS  accuracy  with  the  protection  INS  systems 
have  against  external  noise  [2, 4, 13, 14, 17, 19, 21, 23, 39, 44, 45, 50, 51, 56]. 

Integration  of  navigation  aids  can  further  improve  the  performance  of  the  system.  Due  to 
the  dependence  on  GPS,  it  would  be  useful  to  have  a  system  that  could  determine  the 
level  of  accuracy  and  reliability  of  GPS  transmissions  before  computing  a  navigation 
solution. 
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1.1.4.  Sensor  Integration 

An  extended  Kalman  filter  can  combine  the  position,  velocity,  and  attitude  states 
output  by  the  INS  with  measurements  from  a  GPS  receiver,  barometric  and  radar 
altimeters,  and  a  pseudolite  to  determine  the  best  estimate  of  those  states.  The  Kalman 
filter  in  this  research  uses  tight  integration  to  process  the  measurements.  Tight 
integration  incorporates  raw  measurements  from  all  sensors,  e.g.,  pseudoranges  for  GPS, 
into  a  central  Kalman  filter  as  shown  in  Figure  4.  In  contrast,  loose  integration  provides 
each  sensor  (or  a  selected  subset  of  sensors)  with  its  own  filter  to  process  its 
measurements  before  incorporating  the  measurements  into  the  central  filter  as  shown  by 
Figure  4.  Although  loose  integration  offers  less  complexity  for  each  filter,  tight 
integration  is  preferred  in  this  case,  as  it  requires  the  design  and  tuning  of  only  a  single 
filter  and  has  better  performance  potential  [20]. 

Furthermore,  in  cases  where  fewer  than  four  GPS  signals  are  available  (the  minimum 
to  compute  positions  and  velocities  using  GPS),  the  filter  can  still  function  properly  [27]. 
If  a  separate  filter  is  devoted  to  processing  GPS  measurements,  as  in  loose  integration,  it 
cannot  give  an  unambiguous  estimate  when  there  are  fewer  than  four  measurements 
available.  The  entire  benefit  of  GPS  is  lost  during  those  periods.  When  tight  integration 
is  implemented,  all  measurements  are  delivered  to  a  single  state  estimation  filter.  When 
there  are  fewer  than  four  GPS  measurements  available,  the  filter  continues  to  function, 
deriving  some  advantage  of  having  at  least  partial  information  from  GPS. 
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Best  Estimate  of  T rue  States  Best  Estimate 


Best  Estimate  of  True  States  Best  Estimate 


Figure  5.  Configuration  using  Loose  Integration 
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1.2.  Problem  Definition 


The  core  of  this  work  is  to  examine  the  application  of  the  newly  developed  Modified 
Multiple  Model  Adaptive  Estimation  (M^AE)  algorithm  [33]  to  an  existing  navigational 
system  model  using  a  GPS-aided  INS  in  the  presence  of  external  interference.  This 
algorithm  is  then  applied  to  the  precision  landing  problem  to  determine  its  suitability. 

1.3.  Previous  Research 

A  great  deal  of  research  has  been  conducted  to  develop  an  adequate  precision  landing 
system  (PLS)  using  GPS  [5, 7, 11, 13, 18,  22, 45,  56].  One  significant  effort  is  the 
FAA’s  Local  Area  Augmentation  System  (LAAS),  which  is  based  upon  a  ground-based 
system  broadcasting  DGPS  error  corrections  and  utilizing  airport  pseudolites  [5, 18]. 
Within  AFTT,  there  have  been  two  realms  of  research  applicable  to  PLS.  The  first  relates 
to  the  practical  aspects  of  PLS,  which  include  Johnson  [23]  and  Negast  [38]  who  focused 
on  GPS/INS  navigation  systems.  The  work  of  Gray  [13]  applied  the  system  to  a  PLS  and 
showed  that  ELS  Category  II  could  be  satisfied  using  a  GPS-aided  INS  with  radar 
altimeter  and  pseudolite  measurements  incorporated.  Britton  [7]  follows  this  by  using 
differential  GPS  (DGPS)  measurements  to  show  Category  m  could  be  met.  White  [56] 
investigated  the  performance  a  PLS  in  the  presence  of  GPS  interference  and  spoofing. 
White  implemented  Multiple  Model  Adaptive  Estimation  (MMAE)  [31]  techniques  to 
adjust  the  navigational  system’s  model  to  give  dependable  results.  Advanced 
developments  in  more  theoretical  topics  involving  the  M^AE  architecture  expanded  the 
capabilities  of  MMAE  to  handle  the  task  of  simultaneous  estimation  of  interference  levels 
and  navigation  performance  [33,  55]. 
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1.4.  Scope  of  the  Problem 

The  basis  of  this  work  is  the  expansion  of  the  effort  at  AFIT  towards  investigating 
PLSs  [56].  It  will  primarily  be  the  application  of  the  algorithm  developed  by  Miller 
known  as  Modified  MMAE  (M^AE)  [33]  and  to  extend  the  work  of  White  to  simulate 
varying  GPS  interference  levels  for  a  GPS/INS  PLS  [33].  This  research  uses  a  computer- 
based  simulation,  Multiple  Model  Simulation  for  Optinuil  Filter  Evaluation  (MMSOFE) 
[39],  to  test  the  performance  of  a  model  of  the  PLS.  The  model  of  the  PLS  implemented 
at  AFTT  assumes  the  INS  is  the  primary  navigation  aid. 

Two  configurations  are  investigated  using  various  INS  models  with  two  different  drift 
rates,  combined  with  various  combinations  of  other  navigation  aids  to  determine  the 
effects  on  the  overall  performance  of  the  system.  The  principal  criteria  by  which  to  judge 
the  performance  are  the  FAA  precision  approach  categories.  All  aircraft  navigational 
computers  can  only  use  finite-order  models  of  the  system.  To  reflect  this,  an  elaborate  96- 
state  truth  model  is  used  to  simulate  “real-world”  aircraft  position,  velocity,  and  attitude 
errors,  and  a  reduced-order  13-state  filter  model  serves  as  the  aircraft’s  on-board 
navigation  filter  model  of  sensor  error  characteristics.  The  difference  between  the  states 
of  the  truth  model  and  state  estimates  of  the  filter  model  will  provide  the  error  data  to  be 
weighed  against  the  FAA  categories.  Monte  Carlo  analysis  is  conducted  to  tune  each 
filter  model  for  optimal  performance  and  to  determine  whether  the  configuration  is 
sufficient  for  a  PLS.  Monte  Carlo  analysis  is  conducted  using  MMSOFE,  along  with 
MPLOT.  MPLOT  is  a  software  tool  used  in  conjunction  with  MMSOFE  that  extracts  the 
raw  data  from  the  files  created  by  MMSOFE  to  provide  usable  data  for  plotting  [36]. 
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In  addition,  this  work  serves  as  an  initial  application  of  M^AE  techniques  to  a 
practical  situation.  The  M^AE  concept  has  been  verified  using  truth  and  filter  models  of 
the  same  order  [33].  The  capability  of  M^AE  using  a  more  realistic  and  complete  truth 
model  and  a  reduced-order  filter  model  in  the  simulation  is  accomplished. 

1.5.  Assumptions 

As  in  all  research,  a  number  of  assumptions  must  be  made  to  define  the  limits  of  the 
corresponding  results  and  conclusions.  The  most  crucial  elements  are  the  models  used  to 
depict  the  behavior  of  navigation  equipment  in  the  real  world.  The  96-state  truth  model 
is  assumed  to  represent  how  the  system  would  perform  in  the  real  world,  while  the  filter 
model  serves  as  an  example  of  how  an  operational  navigation  system  would  be 
implemented.  The  simulated  flight  profile  and  GPS  satellite  ephemeris  data  used  are 
generated  by  a  program  called  PROFGEN  [35]  designed  specifically  to  create  the  proper 
files  for  use  with  MMSOFE  [39]. 

The  INS  is  assumed  to  have  a  barometric  altimeter  integrated  with  it,  as  this  is 
commonly  done  to  stabilize  the  INS’s  vertical  readings.  A  feed-forward  configuration  is 
used  with  the  INS,  which  means  the  INS  receives  no  data  to  update  its  current  position. 
Although  a  feedback  configuration  in  which  the  INS  could  accept  position  updates  from 
the  navigation  filter  could  yield  better  accuracy,  if  the  measurements  of  a  sensor  were 
corrupted,  the  INS  would  be  extrapolating  its  states  based  on  an  erroneous  input  which 
can  further  distort  its  results. 

The  presence  of  GPS  interference  is  modeled  as  zero-mean  white  Gaussian  noise 
(WGN)  in  the  measurement  models.  This  is  the  only  interference  to  be  taken  into 
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account;  the  INS  is  considered  to  be  unaffected  by  outside  disturbances,  and  jamming  of 
the  radar  altimeter  is  not  examined.  Pseudolites  can  be  affected  in  the  same  manner  as 
GPS  satellite  signals. 

1.6.  Summary 

This  work  follows  the  progress  of  other  studies  concerning  integrated  GPS/INS 
navigation.  The  intent  is  to  incorporate  the  recently  developed  theory  behind 
simultaneous  parameter  and  state  estimation  with  M^AE  into  this  path  of  study,  PLS 
improvement.  Chapter  2  presents  the  theory  behind  a  critical  tool  of  modem  navigation 
systems,  Kalman  filtering.  It  then  follows  with  a  description  of  parameter  and  state 
estimation  using  MMAE,  and  provides  an  overview  of  the  development  behind  M  AE. 
Chapter  3  presents  the  system  models  implemented  in  the  simulations.  The  structure  of 
the  models  for  each  navigation  system  component  is  described  as  well  as  the  states  and 
measurements  modeled.  Chapter  4  shows  the  results  obtained  through  simulations  and 
the  corresponding  analysis  of  the  performance  observed.  Lastly,  Chapter  5  gives  a 
summary  of  conclusions  based  on  this  work  and  any  recommendation  for  other  paths  of 

study. 
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2.  Theory 


2.1.  Overview 

The  overall  structure  of  the  navigation  system  simulated  in  this  work  is  based  on  a 
modified  MMAE  (M^AE).  Miller  developed  this  method  as  a  solution  to  the  dilemma  in 
simultaneously  estimating  both  the  states  and  the  parameters  of  a  system  model 
accurately  [33].  The  structure  of  M^AE  is  composed  of  two  major  parts:  the  MMAE, 
designed  specifically  for  parameter  estimation,  and  a  single  extended  Kalman  filter 
(EKF)  tuned  precisely  for  state  estimation,  once  provided  the  parameter  estimate  from  the 

MMAE. 

The  MMAE  is  composed  of  a  bank  of  Kalman  filters,  each  one  based  upon  a  different 
value  of  parameters.  Both  the  MMAE  and  state  estimator  receive  measurements  from  the 
external  sensors.  The  MMAE  uses  these  data  to  come  up  with  a  blended  parameter 
estimate,  which  is  passed  on  to  the  state  estimator.  The  state  estimator  then  has  the  best 
model  available  to  render  its  estimate  of  the  states,  based  on  the  incoming  measurements. 
In  this  way,  accurate  state  estimates  can  be  obtained  in  an  environment  where  some  of  the 
parameters  of  the  model  can  change.  The  Kalman  filter  serves  as  the  fundamental 
building  block  of  an  M^AE;  thus,  an  explanation  of  the  M^AE  begins  from  the  most 
elementary  component  proceeding  towards  the  overall  structure. 

This  study  uses  Kalman  filters  to  combine  various  sensor  data  (INS,  GPS,  etc.)  to 
obtain  an  optimal  estimate  of  the  aircraft’s  position  and  velocity.  Kalman  filtering  has 
been  used  considerably  in  navigation  applications  [19, 27]. 
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Figure  6.  M^AE  Structure 


It  is  well  suited  for  this  due  its  ability  to  take  into  account  levels  of  uncertainty  of  states 
and  measurements  within  a  system.  Kalman  filtering  deals  with  systems  that  can  be 
described  by  a  set  of  linear,  time-varying,  stochastic  differential  equations  [24, 30].  In  a 
sampled-data  environment,  such  as  a  computer,  discrete  difference  equations  can  be  used 
with  some  adaptation.  The  stochastic  processes  involved  are  represented  by  white 
Gaussian  noise  (WGN),  usually  zero-mean  unless  a  bias  is  needed.  “White”  indicates 
that  the  value  of  the  process  at  an  instant  of  time  is  independent  of  all  other  values  at  any 
other  time,  i.e.,  “perfect”  randomness.  Gaussian  functions  are  used  since  uncorrelated 
jointly  Gaussian  functions  are  also  independent.  The  strength  values  can  be  adjusted  to 
fit  the  level  of  uncertainty  needed  for  the  models. 

2.2.  Extended  Kalman  Filter 

Originally,  Kalman  filters  were  designed  to  handle  systems  adequately  described  by 
linear  models  [30].  All  models  are  approximations,  and  in  many  cases,  linear  models 
have  deficiencies  that  are  not  negligible.  This  is  often  the  case  when  the  noise  strengths. 
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i.e.,  levels  of  uncertainty,  of  the  model  are  small  values.  Over  time,  the  states  of  the  filter 
can  rely  too  heavily  on  the  system  model  and  become  affected  less  by  external 
measurements.  This  can  result  in  model  behavior  that  strays  away  from  the  actual  real- 
world  situation  [31].  Several  methods  have  been  developed  to  counteract  this.  Some 
restrict  the  variance  of  the  states  from  going  below  a  certain  value  or  have  a  built-in 
minimum  threshold  of  uncertainty.  Others  limit  the  “memory”  of  the  filter,  such  as  in 
Fagin  age-weighting  where  the  effects  of  more  recent  measurements  are  magnified  by 
increasing  the  assumed  noise  variance  of  earlier  measurements  exponentially  backwards 
in  time  [31].  In  cases  where  the  linear  model  cannot  sufficiently  characterize  the  system, 
the  filter  equations  can  be  adapted  to  handle  non-linear  equations.  A  linearized  Kalman 
filter  assumes  that  a  nominal  state  trajectory  (value  of  the  state  vector  over  time)  exists 
and  estimates  the  perturbation  about  the  nominal  value.  The  perturbation  is  defined  as 
the  first-order  term  of  the  Taylor  series  of  the  difference  between  the  state  vector  and  its 
nominal  trajectory.  This  filter  performs  well  unless  the  actual  and  nominal  state 
trajectories  differ  significantly  [31].  To  compensate  for  this,  extended  Kalman  filtering 
(EKF)  relineaiizes  about  the  most  recent  state  update.  Thus,  the  filter  computes  a  new 
nominal  state  trajectory  after  each  update  phase. 

In  the  equations  that  follow,  a  certain  type  of  notation  is  used,  following  a  convention 
adopted  by  Maybeck  [30,  31,  32].  Vector  and  matrix  variables  are  distinguished  from 
scalar  variables  through  a  boldface  type.  Boldface,  lower-case  letters  are  reserved  for 
vectors  and  capital  letters  for  matrices.  Thus,  a,  b,  c  would  be  scalars,  a,  b,  c  would  be 
vectors,  and  A,  B,  C  would  be  matrices.  Variables  are  further  distinguished  by  typeface. 
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Roman  type  indicates  a  deterministic  variable,  e.g.,  a,  b,  c,  while  a  sans-serif  (Helvetica) 
type  indicates  a  stochastic  variable,  e.g.  a,  b,  c. 


2.2.1.  System  Dynamics  Model  Equation 

Consider  the  linear  system  dynamics  model  equation  used  for  Kalman  filtering 

x(0  =  F(0x(0  +  B(0u(0  +  G(t)  w(0  ( 1 ) 


where  x(0  is  the  state  vector,  F(0  is  the  system  dynamics  matrix,  u(0  is  a  vector  of 
deterministic  control  inputs  with  B(t)  as  a  control  input  matrix,  and  w(t)  is  a  vector  of 
zero-mean  WON  with  G(t)  as  a  noise  input  matrix.  Note,  w(t)  is  the  hypothetical  time 
derivative  of  a  Brownian  motion  process  P(0 


W(0  = 


dm 

dt 


(2) 


which  does  not  exist  in  the  real  world  since  w(0  has  an  equal  power  density  over  all 
possible  frequencies  (hence,  the  term  “white”).  However,  if  it  is  assumed  that  a  finite 
band  of  frequencies  (of  concern  to  the  problem)  is  used,  then  w(0  “appears”  as  a  white 
noise  process  within  that  band.  Also,  since  w(t)  is  a  Gaussian  process,  two  statistics, 
mean  and  (auto)covariance  kernel,  are  sufficient  to  completely  describe  how  it  behaves: 

£:{w(0}  =  0  (3) 

E{w(Ow^(t  +  t)}  =  Q(0^(^)  (4) 

Q(t)  is  a  square  matrix  that  can  be  interpreted  as  the  strength  of  the  WGN,  while  6{t)  is 
the  Dirac  delta  function,  also  known  as  the  impulse  function.  With  the  addition  of  w(0, 
X(0  becomes  a  stochastic  variable  itself  with  these  statistics: 

E{x(0}  =  x(0  (5) 
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E[  [x(0  -  x(0][x(0  -  x(Of }  =  P(0 


(6) 


x(?)  is  interpreted  as  the  best  estimate  of  the  state  vector  at  time  t,  and  P(0  is  its 

corresponding  covariance  or  range  of  uncertainty. 

When  the  system  dynamics  cannot  be  expressed  as  a  set  of  linear  equations,  a  non¬ 
linear  function  must  be  set  up  for  the  model: 

x(0  =  f[x(t),u(0,t]  +  G(0w(0  (7) 

where  f[x(r),  u(0,  t]  is  the  non-linear  vector  function  in  terms  of  the  state  and  control 
vectors  and  time.  Also,  in  this  research,  no  deterministic  control  inputs  are  considered, 
thus  B(0  and  u(0  are  zero  and  do  not  appear  from  this  point  on. 

2.2.2.  Measurement  Model  Equation 

The  linear  discrete  measurement  model  equation  used  for  Kalman  filtering  is 

z(r,)  =  H(r,)x(r,)+v(r,)  (8) 

where  z(tf)  is  the  measurement  vector,  H(r,)  is  the  measurement  matrix,  and  v(r,)  is  a 
vector  of  zero-mean,  discrete-time  WGN  representing  the  uncertainty  of  the 
measurements,  or  the  measurement  noise.  It  is  important  to  note  this  equation  is  a 
discrete-time  process  (r,-  instead  of  t)  since  measurements  are  taken  at  specific  instants  of 
time.  Here,  the  measurements  are  considered  to  be  taken  at  regular  time  intervals.  The 
measurement  noise  v(r,)  has  the  following  statistics: 

E{v(r,)}  =  0  (9) 

(10) 
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In  this  description,  R(t,)  is  a  square,  positive  definite  matrix  interpreted  as  the  covariance 
of  the  noise  v(t,).  R(t,)  being  positive  definite  implies  that  all  measurements,  i.e., 
elements  of  V(f,),  are  corrupted  by  noise,  and  no  linear  combination  of  these 
measurements  would  be  free  from  noise.  In  addition,  it  is  assumed  that  the  state  and 
noise  vectors,  x(0,  w(0,  and  v(f,),  are  uncorrelated  with  each  other,  and  being  jointly 
Gaussian  this  means  they  are  also  independent  of  each  other[30]. 

As  above,  when  the  measurement  dynamics  cannot  be  expressed  as  a  set  of  linear 
equations,  a  non-linear  function  must  be  used: 

z(t,)=h[xa,),t,]+v(f,)  (11) 

where  h[x(0,  t]  is  the  non-linear  vector  function  in  terms  of  the  state  vector  and  time. 


2.2.3.  System  Model  Linearization 

Equations  (7)  and  (11)  describe  the  non-linear  model  of  the  system  of  which  a 
Kalman  filter  determines  the  best  estimate  of  the  states  x{t)  and  the  estimation  error 
covariance  P(0.  A  linear  approximation  must  be  made  to  create  a  linearized  Kalman 
filter  to  apply  to  the  model.  A  nominal  state  trajectory  x„(0  V  r  g  7(7  is  the  entire 
period  the  filter  operates)  is  computed  with  initial  conditions  of  xM  =  x„^,  where  the 
noise-free  system  dynamics  equation  is  defined  by 

x„(0  =  f[x„(0.^]  (12) 

where  is  the  same  as  in  Equation  (7).  Nominal,  noise-free  measurements  are  also 
considered  with  the  corresponding  measurement  equation 

z„(/,.)  =  h[x„(t,.),t,]  (13) 

where  is  the  same  as  in  Equation  (11). 
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As  mentioned  in  the  beginning  of  this  section,  a  linearized  Kalman  filter  estimates  a 
perturbation  to  the  nominal  trajectory.  The  state  perturbation  is  found  by  subtracting  the 
nominal  states  from  the  original  states.  This  yields  a  system  dynamics  equation  of 

[x(0-x„a)]  =  {f[x(t),t]-f[x„(0,^]}  +  G(0w(0  (14) 

Now,  a  Taylor  series  expansion  is  performed  on  f[x(r),  t\  about  x„(0  te  T 

=  (X(0-X.(01‘  (15) 

t=0  dx  x(;)=x„(f) 

A  first-order  approximation  may  be  made  by  neglecting  all  the  terms  with  powers  greater 
than  one  (k  >  2) 

f[x(t),t]  f[x„  [x(0-x„(0]  (16) 

x(0=x„(r) 


The  same  procedure  is  applied  to  the  measurement  equation,  with  a  Taylor  series 
expansion  of  h[x(t,),  t]  about  z„(t,)  V  /,•  G  T  yielding 

5z(t,)  =  H[r;x„(t,)]5x(0 + v(0  (19) 

At  this  point,  the  Kalman  filter  equations  can  be  applied,  resulting  in  the  linearized 
Kalman  filter.  It  should  be  noted  that,  unlike  the  basic  Kalman  filter,  the  best  estimates 
of  the  perturbation  states  dk(t)  are  first-order  approximations  and  not  truly  optimal.  In 
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order  to  get  the  estimates  of  the  states  themselves,  the  nominal  state  trajectory  must  be 
added 

x(t )  =  x„  (0  +  dk(t)  (20) 

2.2.4.  Extended  Kalman  Filter  Equations 

In  many  cases,  the  effect  of  neglecting  the  higher  order  terms  of  the  Taylor  series  can 
cause  the  accuracy  of  the  filter’s  predictions  to  stray  over  time.  This  is  caused  when  the 
nominal  state  trajectory  deviates  significantly  from  the  actual  trajectory.  Extended 
Kalman  filtering  solves  this  by  relinearizing  about  the  most  recent  state  estimate  x(t) 
instead  of  the  initial  nominal  state  trajectory  x„it).  Thus,  a  new  nominal  trajectory  is 
created  for  each  cycle  of  the  filter.  Each  cycle  of  a  Kalman  filter  has  two  phases: 
propagation  and  update.  Each  filter  update  occurs  at  regular  intervals  t,  {i  =  0, 1,  2, ...) 
and  of  instantaneous  duration.  The  propagation  phase  takes  place  during  the  intervals 
between  t,-  and  tM.  In  the  EKE  equations  below,  the  following  notation  convention  is 
observed  [30, 31, 32]: 

t\t^  indicates  the  value  of  a  given  variable  at  time  t,  conditioned  on  measurements 
taken  through  time  U  (this  also  represents  the  relinearization  taking  place). 
t:  indicates  the  value  after  the  propagation  phase,  prior  to  the  update  phase 
t*  indicates  the  value  after  the  update  phase,  prior  to  the  next  propagation  phase 

2.2.4.I.  Propagation  Phase 

The  values  of  the  state  estimate  x(t|t,.)  and  state  covariance  P(t|t,.)  are  propagated 
from  t^  to  tT^i  through  the  following  differential  equations: 
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P((|t,)  =  F[t;i((|<,)]F((|(,)+P(tl/,)F’'[nx(r|(,)]+G(()Q(t)G"(() 


where 


F[f;x(f|r,)]  = 


af[x(0,?] 


dx 


with  these  initial  conditions: 


x(tjt,)=x(t;) 

p(dt,)=p(t;) 


(21) 

(22) 

(23) 

(24) 

(25) 


2.2.4.2.  Update  Phase 

After  each  propagation,  the  discrete-time  measurements  z,-  =  z(r,)  are  applied.  The 
updated  values  of  the  state  estimate  x(r,.)  and  state  covariance  P(r,  )  are  computed  through 


the  following  equations: 

K(r,)  =  P(f")H^[tr4(tr)]{H[t,4x(/-)]P(tr)H^[t,;x(t^^ 
x(r;)  =  x(r:)+K(t,){z,  -h[x(r-),/,]} 


where 


H[r,x(r|r,)]  = 


dh[xit),t] 


dx 


x(f)=x(r,-  ) 


and  translating  the  values  after  the  propagation  phase  as: 

x(r-)  =  x(r,.|t,._,) 

P(t:)  =  P(t,.|t,._,) 


(26) 

(27) 

(28) 

(29) 

(30) 

(31) 
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The  values  x(t!)  and  Pit*)  are  then  used  to  start  the  next  propagation  phase  and  to 
recalculate  the  nominal  state  trajectory. 

The  term  {z,.  -h[x(t: ),?,  ]}  in  Equation  (27)  is  known  as  the  measurement  residual, 
and  is  represented  by  r(r,).  The  value  of  r(ti)  represents  the  difference  between  the  actual 
measurements  taken  z/  =  z(t,)  and  the  filter’s  prediction  of  the  measurements 
z(r,)  =  h[x(r,"),t,  ]  based  on  the  best  estimate  of  the  state  vector  prior  to  receiving  z,-.  The 
characteristics  of  r(t,)  show  how  well  the  filter  model  simulates  the  behavior  in  the  real 
world.  In  a  linear  Kalman  filter,  for  an  accurate  model,  r(t,)  will  appear  white  and 
Gaussian  with  a  mean  of  zero  and  a  covariance  of 

A(t,)  =  H[t,.;x(tr)]Par)H^[^,;x(/:)]  +  R(t,)  (32) 

Note  this  term  is  embedded  in  the  expression  for  K(ri)  in  (26).  For  an  extended  Kalman 
filter,  this  description  is  good  only  to  the  first  order. 

2.3.  Multiple  Model  Adaptive  Estimation  (MMAE) 

As  a  Kalman  filter  runs,  the  algorithm  generates  its  prediction  of  the  states  of  the 
system  model,  x(0,  along  with  an  estimation  error  covariance  P(0.  All  of  this  depends 
on  how  well  the  model  itself  depicts  the  system’s  behavior  in  the  real  world.  The 
structure  of  the  model  is  based  on  the  parameters  describing  the  interrelations  of  the 
states.  For  this  work,  these  are  the  system  dynamics  function  f[x(0,  t],  system  noise 
input  matrix  G(t),  measurement  function  h[x(0,  t],  and  the  strength  of  the  system 
dynamics  noise  Q(0  and  covariance  of  the  measurement  noise  R(t,)-  The  filter  can 
encounter  difficulties  when  any  of  these  quantities  change  from  values  assumed  by  the 
filter.  Thus,  it  would  be  beneficial  to  have  a  filter  that  could  adapt  to  such  changes.  A 
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technique  known  as  multiple  model  adaptive  estimation  (MMAE)  handles  changes  in 
system  characteristics  through  a  parallel  bank  of  K  individual  Kalman  filters  (where  K  is 
the  number  of  possible  parameter  values  ajt  taken  into  account),  each  calculating  its  own 
best  estimate  of  the  state  vector  with  P^(t.)  (k  =  1, 2, K).  When  the 
measurements  z,-  are  given  to  the  filters,  they  generate  a  set  of  residuals  Tidti).  These 
residuals  are  processed  by  a  hypothesis  conditional  probability  computation  algorithm  to 
quantify  how  well  each  filter  matches  the  real  world  circumstances.  This  is  expressed  as 
a  conditional  probability,  pkitd,  by  the  recursive  equation 

^  /z«,  )la,Z(o., )  1®  J  ’  ^i-l )  Pj  ^^i-1  ^ 

i=i 

where  a*  is  the  value  of  parameters  associated  with  the  k'^  filter,  Z/.i  is  the  realization  of 
the  variable  of  Z(t,.i),  which  is  in  turn  a  vector  of  cumulative  measurements  from  the 


initial  time  to  to  ti-i: 


and  /j(,  )|a  z(,.  ,)(Z/|a*  is  a  Gaussian  density  function  with  a  mean  of  rt(t,)  and  a 


covariance  of  Ak(ti): 


where 


(2;r)^|A^(t,)l' 
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and 


(37) 


The  quantity  rji(t,)  is  an  m-dimension  vector,  i.e.,  m  scalar  measurements  are  taken  at 
each  measurement  update  interval.  These  conditional  probabilities  are  then  multiplied  by 
their  corresponding  state  estimates,  which  are  summed  to  produce  a  blended  estimate  of 
the  state  vector  x(r,.).  Figure  7  shows  a  block  diagram  of  the  MMAE  [31]. 

The  conditional  probability  is  best  interpreted  as  pkiti)  =  Prob[a  =  a*  |  Z(ti)  =  Z,],  i.e., 
the  probability  that  the  parameter  value  ajt  is  the  correct  one  for  the  model,  based  on  the 
measurements  collected  up  to  this  point  Z(t,)  =  Z,-. 


26 


The  blended  estimate  is  then  the  weighted  mean  of  =  1, 2, K\ 

x(r,.)  =  £;{x(t,)|Z(t,)  =  Z,}  =  (38) 

The  corresponding  covariance  P(^/)  of  the  blended  solution  is 

P(^,)  =  2{P*(f,)+[x*(^/)-x(t,)][x,(/,)-x(t,)f  }A(t,)  (39) 

k^\ 

The  underlying  concept  of  the  MMAE  is  that  values  of  the  residuals  determine  which 
filter  best  represents  the  real-world  system  at  the  current  time,  i.e.,  the  best  filter  has  the 
smallest  scaled  residuals.  These  residuals  are  scaled  by  the  inverse  of  the  filter-computed 
residual  covariance,  as  seen  in  Equation  (37).  The  corresponding  probability  of  that  filter 
should  increase  over  time  and  the  other  filters’  probabilities  should  decrease,  assuming 
the  system  environment  remains  constant.  Prior  to  using  the  MMAE,  each  filter  should 
be  tuned  (values  of  the  dynamics  pseudonoise  strength  matrix  Q  adjusted)  for  its  best 
performance  by  running  the  algorithms  under  the  conditions  of  each  of  K  parameter 
values  being  the  best  or  “true”  values  in  the  real  world.  Thus,  the  filter  would  be 
tuned  for  best  performance  when  the  “true”  parameters  are  a*. 

An  inherent  trade-off  problem  exists  in  an  MMAE  between  state  estimation  (x  )  and 
parameter  estimation  (a).  The  problem  is  that  the  only  way  an  MMAE  can  ascertain 
what  the  system  parameters  should  be  is  through  the  measurements  it  receives.  When  an 
MMAE  is  designed  to  yield  good  state  estimates  (i.e.,  accurately  predicting  the  system’s 
actual  states  X(t,)  given  all  previous  actual  measurements  Zm  over  a  period  of  time)  each 
filter  is  tuned  individually  to  provide  accurate  state  estimates  for  its  own  parameter  value 
a<:.  If  the  filters  are  tuned  conservatively,  the  residuals,  (t, )  =  z,  -  h[X;^  (t,. ),  t,-  ] ,  can 
appear  similar  to  each  other.  If  all  the  residuals  tend  to  have  nearly  the  same  magnitude. 
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then  the  MMAE  favors  (assigns  the  highest  probability  to)  the  filter  with  the  least  value 
of  |A,|,  as  shown  by  Equations  (33)  and  (35)  -  (37).  Although  A^  is  used  to  represent  the 

covariance  of  the  residuals  (assumed  Gaussian),  the  value  of  [A^^l  is  independent  of  the 
residuals  themselves  and  how  well  the  ^***  filter  fits  the  parameter  environment.  The 
conditional  probability  of  each  filter  pkit.)  is  dependent  upon  its  residuals  through  Lkitd 
from  Equation  (37).  If  the  scaled  values  of  the  residuals,  Tkiti),  are  similar,  then  Li(t,)  - 
L2(ti)  =  . . .  «  L/f(t/)-  This  means  that  pkiti)  is  how  highly  dependent  on  Pkitd  (2-32)  which 
has  as  its  only  variable,  |A;t(t,)| .  A  strong  dependence  on  |At(t,)|  can  distort  the 
blending  process  since  there  will  be  one  filter y  e  {1,2,  such  that  |A^(t,)|  < 

|A^  (f.)|  \f  k^j.  If  the  parameter  environment  does  not  change  for  a  period,  pjiti)  will 

approach  unity  while  all  other  pk{ti)  will  approach  zero.  As  this  happens,  one  filter  out  of 
the  bank  becomes  favored  by  having  most  of  the  probability  weighting.  This  may  be  the 
“best”  filter  out  of  the  choices  available,  though  it  may  not  have  the  best  parameters  out 
of  the  entire  parameter  space. 

However,  the  performance  of  the  algorithm  relies  upon  significant  differences 
between  the  residual  characteristics  of  each  elemental  filter.  Because  of  this,  it  is 
important  to  avoid  adding  too  high  levels  of  dynamics  pseudonoise  Q  during  filter  tuning. 
Although  a  conservative  tuning  philosophy  is  used  to  keep  filter  estimates  from  diverging 
considerably  from  the  truth,  conservative  tuning  tends  to  mask  the  differences  between 
good  and  bad  models  [28,  33]. 

Another  concern  occurs  when  one  of  the  elemental  filters  achieves  a  conditional 
probability  of  zero.  Whenever  this  occurs,  that  filter’s  probability  remains  zero 
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indefinitely  due  to  the  recursive  nature  of  Equation  (33).  Thus,  the  filter  is  eliminated 
from  the  possible  choiees  considered  by  the  MMAE  in  the  future.  Even  if  the  parameter 
conditions  change  such  that  the  elemental  filter  that  had  zero-probability  now  represents 
the  best  model  of  the  real-world  system,  that  filter  cannot  have  any  probability  weight 
assigned  to  it.  One  solution  to  prevent  this  from  happening  is  to  establish  an  artifieial 
minimum  threshold  value  that  pk  can  attain,  i.e.,  min(pii:)  [1, 31].  Consequently,  this 
means  no  elemental  filter  can  have  =  1  even  if  it  is  the  exact  match,  though  this  should 
not  produce  any  significant  problems  in  determining  a  solution  [31]. 

The  advantage  MMAE  provides  to  Kalman  filtering  is  that  it  handles  situations  where 
a  single  system  model  cannot  sufficiently  simulate  the  behavior  of  an  actual  system 
confronted  with  varying  parameters.  In  addition  to  the  states,  the  MMAE  algorithm  can 
be  used  to  determine  the  best  estimate  of  the  parameters  themselves  a(t,  )  based  on  all  the 
measurements  collected  Z(t,)  =  Z,.  The  blended  estimate  of  the  parameter  set  is 
calculated  in  the  same  manner  as  with  the  state  vector: 

a{(,)  =  £{a(0|Z((,)  =  Z,}  =  2a, ((,)?,«,)  (40) 


The  conditional  covariance  of  a(t,)  of  the  blended  parameter  estimate  is 


P,(r.)  =  £{[a(t,)-a(0][a(t,)-a(t,)r|Z(t,)  =  Z,} 
=  S[a,  ih)  -a(t,)][a*  (?,)  -a(t,)f  p*  (t,) 


(41) 


When  the  filters  aio  tuned  for  parameter  estimation,  it  is  done  such  that  the  residuals 
tend  to  have  values  sufficiently  distant  from  each  other.  This  helps  to  create  a  situation 
such  that  if  the  actual  parameters  do  not  exactly  or  closely  match  those  of  one  of  the 
filters,  a  blending  of  several  filters  will  provide  an  accurate  estimate  of  the  parameters. 
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However,  the  method  to  make  the  residuals  more  distinct  from  each  other  entails 
applying  smaller  values  of  the  Kalman  filter  gains  to  reduce  the  influence  of  the 
measurements  z,-.  When  the  measurements  have  less  influence,  the  filter  tends  to  rely  too 
heavily  on  the  propagated  state  estimates  x(tr) .  This  means  the  actual  state  x(f,)  can  be 
misrepresented  over  time.  Therefore,  new  techniques  have  been  developed  to  give 
accurate  estimates  of  both  the  states  and  the  parameters  simultaneously  [33]. 

2.4.  Modified  MMAE  (M^AE) 

The  M^AE  architecture  combines  MMAE  and  Kalman  filter  techniques  for 
simultaneous  state  and  parameter  estimation.  Under  this  architecture,  an  MMAE  serves 
as  the  parameter  estimator  and  is  designed  and  tuned  for  estimating  the  system  s 
uncertain  parameters  accurately.  It  is  optimized  for  distinguishing  among  several 
possible  hypothesized  operating  conditions  dictating  the  parameters  of  the  system.  The 
separate  single  state  estimator  within  the  M^AE  algorithm  is  designed  and  tuned  to 
provide  accurate  state  estimation,  conditioned  on  the  measurements  and  knowledge  of  the 
parameters  provided  by  the  parameter  estimator. 

Several  assumptions  regarding  the  system  must  be  accepted  before  an  M  AE  can  be 
applied  to  the  model.  Most  importantly,  the  parameters  to  be  estimated  lie  within  a  finite 
predefined  parameter  space.  The  elemental  filters  of  the  MMAE  are  then  based  on  a 
subset  of  discrete  parameter  values  chosen  from  the  parameter  space  [31].  In  addition, 
the  variable  parameters  change  more  slowly  than  the  system  s  states,  thus  estimates  of 
these  parameters  can  make  use  of  any  prior  information  available  as  to  how  they  vary 
over  time.  As  for  designing  an  M^AE  architecture.  Miller  describes  a  straightforward 
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method  to  analyze,  tune,  and  predict  the  system’s  performance  before  conducting  a  full- 
scale  Monte  Carlo  analysis.  Figure  8  shows  a  flowchart  of  the  performance  analysis  tool. 

First,  a  discretized  parameter  set  is  established  using  an  algorithm  designed  by 
Sheldon  to  determine  the  parameters  of  each  elemental  filter  of  the  MMAE  [48, 33].  If 
appropriate,  inter-residual  distance  feedback  (IRDF)  techniques  developed  by  Lund  [28] 
are  applied  to  make  each  elemental  filter  appear  more  distinct  from  one  another.  Next, 
the  state  estimation  Kalman  filter  is  designed  using  techniques  described  by  Maybeck 
[30, 31].  It  is  important  to  note  that  the  MMAE  will  provide  the  values  of  the  variable 
parameters  to  the  state  estimator  based  on  the  incoming  measurements.  Then,  the 
MMAE  and  state  estimator  are  coded  in  a  software  simulation.  A  single  Monte  Carlo  run 


Figure  8.  MMAE  Structure 
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is  performed  on  the  MMAE  to  generate  the  values  of  the  parameter  estimate,  its 
covariance  determined  by  the  elemental  filters,  and  the  probabilities  associated 
with  each  elemental  filter  pkiti).  Finally,  an  approximate  covariance  analysis  is  conducted 
to  verify  if  the  design’s  performance  meets  the  desired  specifications.  If  so,  a  thorough 
Monte  Carlo  analysis  is  conducted  on  the  M^AE.  If  not,  the  MMAE  and  the  state 
estimator  are  modified  iteratively  to  solve  the  discrepancies. 


2.4.1.  Parameter  Space  Discretization 

The  M^AE’s  MMAE  is  designed  for  accurate  parameter  estimation.  When  designing 
the  MMAE  for  parameter  estimation,  a  parameter  space  must  be  defined.  The  parameter 
space  is  the  range  of  values  the  uncertain  parameter  to  be  estimated  can  assume.  Once 
this  is  determined,  the  next  issue  is  the  placement  of  the  elemental  filters  to  span  the 
space.  In  the  past,  several  ad  hoc  methods  were  used  to  determine  the  placement  or 
discretization  of  the  parameter  space,  e.g.,  equal  spacing,  exponential  spacing,  etc.  [15, 
16,  31, 46, 49].  Sheldon[48, 49]  thought  a  more  systematic  design  approach  was  required, 
thus  he  developed  a  method  to  assign  the  placement  of  the  filters  optimally.  His  method 
involves  minimizing  a  cost  function,  C,  which  is  the  average  of  the  mean-squared 
estimation  error  taken  over  the  parameter  space: 


j  £;{[a(0-a(0f  W[a(0-a(0]}^^a 


(42) 


where 

^da=  j  ...  j  ^da^  da^ ...  da^  (43) 
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(44) 


a, 

^2 

a=  . 

and  W  is  a  weighting  matrix  chosen  by  the  designer  to  emphasize  specific  states,  N  is  the 
number  of  scalar  parameters  (dimension  of  a),  and  X  is  the  bounded  region  of  the 
parameter  space  within  which  the  parameters  can  take  on  values. 

The  parameter  space  discretization  and  placement  of  the  elemental  filters  is 
accomplished  before  implementing  the  MMAE.  This  entails  deciding  which  parameters 
can  change,  what  the  appropriate  parameter  spaces  are,  and  where  the  elemental  filters 
should  be  placed  within  the  parameter  space.  Sheldon  developed  a  five-step  algorithm  to 
minimize  a  cost  function  over  a  parameter  space  numerically.  The  first  three  steps 
consist  of  constructing  the  truth  and  filter  models  to  represent  the  system,  deciding  the 
number  of  filters  to  be  used,  and  determining  the  cost  function  (parameter,  state,  or 
control,  but  in  the  M^AE,  the  MMAE  is  only  used  for  parameter  estimation)  of  the 
parameter  set  to  be  employed.  The  fourth  step  is  the  core  of  the  process.  The  basic 
purpose  is  to  apply  a  numerical  integration  technique  to  evaluate  the  value  of  the  cost 
function,  C,  over  the  parameter  space  S .  Assuming  the  parameter  set  remains  constant 
for  a  given  problem,  only  the  numerator  of  Equation  (42)  needs  to  be  evaluated,  and  with 
W  often  chosen  as  a  diagonal  matrix,  it  can  be  expressed  as: 

j  E{[a{t)  -  a{t)f  W[a(0  -  a(0]}^^a  =  j  tr( W£:{[a(0  -  a(r)]  [a(r)  -  })da  (45) 

where  tr(»)  is  the  trace  of  a  square  matrix  (sum  of  its  diagonal  elements).  For  numerical 
integration,  K  is  divided  into  a  number  of  discrete  intervals.  At  one  point  of  each 
integration  interval,  a  value  for  a(t,  )  is  calculated  in  the  following  manner: 
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Sheldon  defines  a  transformation  matrix  T  that  converts  true  states  (indicated  by  the 
subscript  T)  x-riti)  to  filter  states  (no  subscript  T)  x(t,)  in  this  way: 

x(f,)  =  TXT(t,)  (46) 

where  x(t,)  is  defined  by  the  discrete-time  system  equation: 

x(t,,,)  =  «D(t,,„t,)x(t,)  +  G,(t,)w,(t,)  (47) 


and 

Ga  (r,  )w,  (t, )  =  j'"‘  >  h  )G(T)i^P(T)  (49) 

and  Wd(t,)  is  zero-mean  WGN  with  strength  (covariance)  of  Qd(^i): 

Q,«,)  =  (50) 

The  following  equation  is  then  solved  iteratively,  where  the  subscript,  k,  denotes  the 
it**'  filter  corresponding  to  ajt  (^=  1,  2, ...,  /iT)  (T,  as  a  subscript,  still  denotes  the  truth 
model): 

r*(n  +  l)  =  Yr*(n)Y^+GoQoGj  (51) 


where 


Y  = 


i^,(I-K,H,)  (0,T-^T^>t)-<I>,K,(HJ-Ht) 
0  <i>T 


Go  - 


-TG„ 


'dX 


0 


Qo  = 


Qdx  0 
0  Rt 


(52) 

(53) 

(54) 
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The  iteration  ends  when  Tkin)  approaches  a  steady-state  value,  i.e.,  limr^(n) ,  (or 

practically  when  Tk{n+1)  ~  Tkin)),  and  is  represented  as  F” .  The  proximity  of  the 
filter  is  defined  as  [49]: 

4=ta|A.|+tr{Ar'[[H.  H.T-H,]rr[H,  H.T-hJ+R,]}  (55) 

If  the  MMAE  has  sufficient  conditions  to  converge  in  the  Baram-sense  [3],  it  will 
converge  toward  the  f"  filter  governed  by: 

Ij  =min{^Jfc  =  l,  2,  ....  K}  (56) 


Then  the  parameter  is  given  by 

^iti)  =  'Z^kPk(h) 

k=l 

=  {[l-/i:p^Ja//,)}  +  Xa.a;K 


(57) 


Jt=l 


where  pimn  is  the  pre-defined  minimum  conditional  probability  a  filter  can  attain,  and  aj  is 
the  value  of  a(r/)  according  to  they*  filter. 

The  fifth  step  is  to  use  the  values  of  a(t,.)  computed  above  in  a  vector  minimization 

of  C(a).  This  yields  the  optimal  values  of  for  each  filter. 

This  procedure  applies  to  steady-state,  constant-gain  (i.e.,  K(r,)  does  not  vary  over 
time)  Kalman  filter.  This  is  adequate  for  time-invariant  systems  with  stationary  noise.  In 
cases  involving  astable  or  unstable  system  model,  the  Kalman  filters  of  the  MMAE  do 
not  achieve  a  steady-state  condition.  This  motivates  the  use  of  a  finite  horizon,  the 
selecting  of  a  period  where  the  system  parameters  can  be  considered  time-invariant  [33]. 
This  research  uses  an  INS  with  time-varying  parameters  in  its  system  dynamics  matrix 
and  an  inherent  instability  in  estimating  vertical  positions  (a  common  trait  of  INSs). 


35 


Thus,  a  finite  horizon  is  implemented  in  the  software  in  the  iterative  solution  of  the 
Equation  (51). 


2.4.2.  Inter-Residual  Distance  Feedback  (IRDF) 

Another  concern  in  parameter  estimation  is  to  keep  the  values  of  Lt(r,)  sufficiently 

distinct  to  avoid  overdependence  on  |Aj^(r,)| .  To  accomplish  this,  Lund  [28, 29,  33] 
defines  a  quadratic  term  JkjiU)  as  the  squared  distance  measure  of  the  difference  between 
residual  TkiU)  and  residual  rj(r,): 

=  <58) 

where  rt/r,)  =  r^Cr,)  -  r/r,),  j  ^  k,  is  the  inter-residual  distance  between  the  A:*  filter  and 
the/**  filter,  and  Wkj  is  a  positive  definite  scaling  matrix.  The  plan  is  to  maintain  Jkjitd 
above  a  specific  limit  value  7o(r,)  >  0  by  adjusting  the  filter  gains.  A  common  method  to 
accomplish  this  involves  adjusting  the  filter  gains  by  varying  the  dynamics  noise 
strengths  downward.  In  the  continuous-time  method  developed  by  Lund,  Qk  is 
replaced  by 

Q;  (0  =  r]it)Qk  (t),  Vit)  e  [ry^in  4]  (59) 

Assuming  Qk  represents  the  dynamics  noise  for  the  filter  tuned  for  best  estimation 
when  its  assumed  parameter  value  is  true,  the  modulating  parameter  r]{t)  has  an  upper 
limit  of  1,  since  ry(r)  >  1  would  decrease  the  distinguishability  of  the  residuals  [28, 29, 
33].  A  lower  limit  rjmin  ^  0  is  set  to  maintain  the  stability  of  the  filter.  In  addition,  Lund 
defines  the  continuous-time  derivative  of  r]{t)  as: 


^  I ^Ukjit) -  4o(0]>  Cond.  I 
dt  ~  [0,  Cond.  2 


(60) 
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where 


Cond.  1:  i](t)  E[ri^„,l] 

Cond.2:  =  AND  J,{t)]<0  OR  (61) 

77(0=  1  ANDa4(0-/o(0]>0 

The  constant  ^  provides  proper  attenuation  of  the  noise  modified  by  rjit).  Ad  hoc 
methods  are  used  to  determine  the  values  of  77inin,  and  Jo-  Lund  suggests  choosing  ^ 
such  that  l/(^  is  greater  than  the  largest  time  constant  of  the  elemental  filters,  and  77nun  =  0 
if  the  system  maintains  its  stability  [28, 29].  For  Jo,  one  ad  hoc  method  involves 
performing  a  sample  run  before  applying  the  IRDF  techniques  [29,  33].  An  initial  value 
for  Jo  is  chosen  based  on  the  actual  values  of  Jkj  observed,  e.g.,  the  mean  value  of  Jkj. 
These  values  are  then  “tuned”  in  further  simulations  to  verify  the  system  performs  within 
its  requirements.  However,  there  is  a  shortcoming  in  this  method  in  that  the  modulated 
system  noise  strength  (t) ,  being  a  function  of  time,  is  computed  in  real-time  and 
cannot  be  determined  in  advance. 

As  a  way  around  this,  Lund  recommends  modulating  the  quantity  Kk(t)Tic(t)  instead  of 
Qk(t)  to  simplify  this  for  linear  models: 

K;(0r,(0  =  t?(0K*(0r*(t),  t](Oe[77^„,l]  (62) 

The  primary  benefit  of  this  is  that  the  Kalman  filter  gains  K^(0  can  be  precomputed 
and  only  r]{f)  is  computed  in  real-time  [29, 33].  In  addition,  modulating  Kjt(r)  yields 
quicker  adaptation  responses  than  with  Qk{t),  which  relies  on  the  filter  state  covariance 
Pjfc(0  equations  transient  effects  to  subside  before  the  filter  gains  can  be  changed  [29, 33]. 
Conversely,  Lund  states  this  approach  has  less  benefit  for  extended  Kalman  filter  and 
other  higher  order  filters  where  their  filter  gains  are  computed  in  the  process.  In  this 
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research  effort,  the  MMAE  filters  implemented  have  reasonably  distant  values  when  the 
parameter  space  is  discretized  by  means  of  Sheldon’s  algorithm.  Also,  many  of  the  terms 
of  the  system  dynamics  noise  matrix  Q  for  the  truth  model  and  tuned  filter  model  are  very 
small  or  zero.  For  these  reasons,  IRDF  was  not  pursued  in  this  work,  although  it  may  be 
of  significant  interest  to  other  potential  research  areas  [29,  33]. 

2.5.  Summary 

The  theory  and  approaches  described  in  this  chapter  are  intended  to  provide  insight 
into  the  methods  used  to  implement  system  models  for  test  simulations.  Due  to  the  extent 
Kalman  filtering  is  applied  to  navigation  problems,  a  good  understanding  of  its  workings 
can  be  invaluable  in  designing  such  systems.  One  of  the  motivations  for  this  work  is  to 
apply  the  newly  developed  M^AE  techniques  [33]  to  an  actual  problem,  thus  improving 
the  capabilities  of  parameter  estimation  (prediction  of  the  system’s  working 
environment),  and  state  estimation  (prediction  of  the  values  of  concern,  such  as  altitude). 
In  this  work,  state  estimation  in  the  face  of  changing  parameters  for  the  precision  landing 
problem  is  the  core  matter. 
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3.  Methodology 


3.1.  Overview 

The  fundamental  elements  of  this  type  of  research  effort  are  the  application  and 
implementation  of  accurate  system  models  and  simulations  for  the  problem  at  hand.  This 
chapter  presents  a  comprehensive  description  of  the  truth  and  filter  models,  integration 
methods,  and  simulation  techniques  employed.  An  overall  description  of  the  system  is 
provided  at  the  start.  Next,  the  truth  and  filter  models  are  developed  for  each  navigation 
subsystem,  including  a  complete  description  of  their  roles.  Finally,  the  approaches  to 
perform  simulations  of  the  problem  are  outlined.  The  next  chapter  then  gives  a 
breakdown  of  the  results  of  the  simulations. 

3.2.  System  Description 

The  navigation  system  to  be  examined  consists  of  several  elements.  The  basis  for  the 
system  is  the  M^AE.  Its  role  here  is  to  provide  the  best  estimates  of  the  error  states  of 
the  INS  rather  than  the  true  states.  True  states  are  the  actual  values  of  position,  velocity, 
attitude,  etc.  This  should  not  be  confused  with  the  truth  model,  which  is  a  system  model 
that  is  assumed  to  depict  the  actual  behavior  of  the  states,  true  or  error.  Error  states  are 
distinguished  from  true  states  in  that  they  are  the  correction  values  of  position,  velocity, 
etc.  that  should  be  subtracted  from  the  states  output  from  the  INS  to  provide  the  best 
navigation  estimates.  This  is  necessary,  since  the  INS,  in  the  configuration  being 
considered,  does  not  have  its  states  reinitialized  to  reflect  an  update  on  the  actual  position 
of  the  aircraft  according  to  the  measurement  devices.  Before  the  entire  system  is  first 
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activated,  the  initial  states  of  the  INS  are  carefully  calibrated  to  reflect  the  craft’s  current 
position,  and  the  M^AE  error  states  are  assumed  to  be  zero.  After  activation,  the  state 
estimator  predicts  how  far  INS  readings  deviate  from  the  actual  (true)  position,  attitude, 
and  velocity  of  the  craft,  hence  the  term  error  states. 

To  compute  its  best  estimate  of  these  error  states,  the  M^AE  receives  information  from 
several  sources.  These  consist  of  the  INS  with  altitude  estimation  {vertical  channel) 
aiding  from  a  barometric  altimeter,  pseudoranges  from  four  GPS  satellites  and  one 
pseudolite,  and  relative  altitude  from  a  radar  altimeter.  The  GPS  signals  and  altimeter 
readings  are  considered  as  measurement  values  for  the  M^AE.  The  truth  model  uses  62 
error  states  to  describe  how  these  systems  function  in  the  real  world  [7, 13,  38, 56]. 

Then,  the  reduced-order  filter  model’s  performance  is  compared  to  the  truth  model’s 
performance.  The  filter  model  tracks  13  error  states,  representing  the  navigation  system 
on  board  the  aircraft.  Figure  9  shows  a  layout  of  the  relationships  between  subsystems. 
For  the  purposes  of  simulation,  a  flight  profile  trajectory  was  generated  in  advance  by  the 
program  PROFGEN.  The  GPS  satellite  data  used  are  from  a  file  recorded  4  May  1991 
[13]. 

The  simulations  compare  the  performance  of  the  13-state  reduced-order  filter  to  the 
62-state  truth  model.  The  (error)  states  of  the  filter  model  are  a  proper  subset  of  those  of 
the  truth  model,  and  were  chosen  to  be  those  having  the  greatest  influence  over  the 
desired  navigation  data  output  of  the  complete  system  [38].  The  truth  model  makes  use 
of  39  states  representing  the  INS  and  23  states  for  the  GPS  data.  From  this  group  of 
states,  11  INS  states  and  2  GPS  states  were  selected  to  be  implemented  into  the  filter 
model. 
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Figure  9.  Navigation  System  Block  Diagram 


3.3.  Inertial  Navigation  System  Models 

The  INS  model  used  in  this  work  is  based  on  the  Litton  LN-93  strapdown,  wander- 
azimuth  INS  using  accelerometers  and  ring-laser  gyroscopes  (RLG)  to  detect  the  changes 
in  the  motion  and  orientation  of  the  aircraft.  Readings  from  these  instruments  are 
translated  into  the  various  states  of  the  model  including  position,  velocity,  and  attitude. 
Litton  developed  a  sophisticated  error  state  model  using  93  states  to  describe  the 
characteristics  of  the  LN-93  INS  [26].  This  error  state  model  provides  an  accurate 
depiction  of  the  actual  performance  characteristics  of  the  INS  and  has  been  used  in 
several  research  projects  at  the  Air  Force  Institute  of  Technology  (AFTT)  [7, 13,  38,  56, 
57].  The  error  state  model  has  been  implemented  into  the  MMSOFE  program  [39,  36]  to 
conduct  the  Monte  Carlo  simulations. 
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The  error  states  of  the  LN-93  Sx.  can  be  expressed  as  a  93-state  vector  that  can  be 
separated  into  six  categories  or  “subvectors”. 

5x  =  [5xJ^  <5X2  5X3  (5xJ  5X5  5Xg]^ 

5xi  has  13  elements  representing  the  commonly  used  nine  Pinson  error 
states:  position,  velocity,  and  attitude  in  three  dimensions;  and  four 
errors  in  the  vertical  channel  of  the  INS  [43]. 

6X2  has  16  elements  representing  the  exponentially  time-correlated  errors  of 
the  gyros,  accelerometers,  and  barometric  altimeter,  and  “trend”  states, 
all  of  which  are  treated  as  first-order  Gauss-Markov  processes. 

5X3  has  18  elements  representing  gyroscopic  bias  errors,  which  are  treated  as 
random  constants. 

5X4  has  22  elements  representing  accelerometer  bias  errors,  which  are 
treated  as  random  constants. 

5X5  has  six  elements  representing  accelerometer  and  initial  thermal 
transients,  all  of  which  are  treated  as  first-order  Gauss-Markov 
processes. 

5x6  has  18  elements  representing  gyroscopic  compliance  errors,  treated  as 
biases. 


(63) 
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The  differential  system  equation  for  the  E^S  truth  model  is 


5x  =  F5x  +  Gw 
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A  complete  description  of  the  submatrices  within  F  is  given  in  Appendix  B. 
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3.3.1.  INS  Dynamics  Truth  Model 

Although  this  93-state  model  provides  the  most  accurate  depiction  of  the  behavior  of 
the  LN-93  INS,  there  is  a  considerable  computational  burden  involved  by  implementing 
this  model  directly.  The  work  of  Negast  provides  a  reduced-order  version  of  this  and 
demonstrates  it  retains  sufficient  fidelity  to  the  original  for  the  types  of  simulations 
conducted  here  [38].  The  reduced-order  model,  which  is  incorporated  into  this  work,  is: 
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The  states  of  this  model  are  a  validated  subset  of  93-state  vector,  though  they  are  not 
the  first  39  states  of  that  vector.  Thus,  the  submatrices  within  this  equation  are  different 
from  those  in  the  previous  equation  since  the  error  state  subvectors  do  not  correspond 
exactly  to  those  above.  The  relationship  between  the  states  of  the  two  models  is  provided 
in  Appendix  B. 
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3.3.2.  INS  Dynamics  Filter  Model 

The  model  representing  a  typical  implementation  in  a  navigation  computer  installed 
on  an  aircraft  would  tend  to  be  less  complex  than  the  full-scale  truth  model.  Processing 
full-order  models  requires  extensive  computing  capacity.  Therefore,  reduced-order 
models  economize  the  limited  computation  resources  normally  available  on  board  an 
aircraft  and  provide  quick  processing  of  results  in  real  time.  This  encourages  the 
development  of  a  further  reduced-order  filter  model.  The  filter  model  for  the  INS  used  in 
this  research  consists  of  1 1  states:  three  position  alignment  errors,  three  velocity  biases, 
three  attitude  alignment  errors  (related  to  the  nine  Pinson  error  states)  [43],  and  two  states 
for  vertical  channel  stabilization.  Since  the  ultimate  goal  is  to  obtain  the  most  current  and 
accurate  information  about  the  position  and  velocity  of  the  craft,  the  states  selected  are 
those  having  the  greatest  influence  on  these  quantities.  All  of  these  states  are  elements  of 
5xi.  Thus,  the  elements  of  the  system  dynamics  matrix  form  a  subset  of  the  elements  of 
Fii.  To  compensate  for  the  omission  of  states,  the  values  of  the  strength  of  the  dynamics 
noise  matrix  Q(0,  where  £:{w(0w'^(r  +  t)}  =  Q(/)<5(t),  t=  0,  are  tuned,  i.e.,  adjusted,  so 
that  the  behavior  of  the  filter  equation  closely  resembles  that  of  the  truth  model. 

Appendix  B  gives  a  further  description  of  these  states  and  Appendix  C  lists  the  tuning 
values. 

3.3.3.  INS  Measurement  Filter  Model 

The  measurement  device  directly  incorporated  into  the  INS  is  the  barometric 
altimeter.  The  barometric  altimeter  provides  signals  to  counteract  the  inherent  instability 
of  the  vertical  channel,  i.e.  altitude  determination,  of  the  INS.  Since  the  M^AE  uses  error 
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states,  the  value  of  the  true  altitude  (unknown)  must  be  eliminated  by  taking  the 
difference  between  the  altitude  predicted  by  the  INS  hm  and  the  reading  of  the 
barometric  altimeter  hsar-  The  value  of  INS-predicted  altitude  hjNs  is  considered  to  be  the 
sum  of  the  true  altitude  hj  and  the  INS  error  in  altitude  above  the  reference  ellipsoid  of 

j 

the  earth  (WGS-84  [9]),  Bi.  The  reference  ellipsoid  is  used  since  barometric  altimeters 
are  calibrated  to  measure  altitude  above  mean  sea  level,  and  not  the  actual  terrain.  The 
barometric  altimeter  reading  is  viewed  as  the  sum  of  the  true  altitude  h-j,  the  time- 
correlated  error  in  the  altitude  reading  dhs,  and  a  random  noise  in  the  measurement  Vgan 
an  element  of  the  measurement  noise  vector  V.  By  taking  the  difference  of  the  two 
altitude  readings,  the  value  of  the  true  altitude  is  eliminated,  thus: 

^Bar  ~  ^INS  ~  ^Bar 

=  [lh+Sh]-[h,+Sh,-v,J  (66) 

=  Sh-Shs+Vg^^ 

The  random  noise  Vsar  is  considered  to  be  zero-mean,  discrete-time  WGN,  thus  it  can 
take  on  a  positive  or  negative  sign.  The  sign  of  VBar  was  chosen  to  yield  a  positive  one  in 
the  final  form  of  the  equation  to  emphasize  the  addition  of  measurement  noise.  Two 
states  from  the  filter  model  appear  in  this  equation:  the  INS  error  in  altitude  above  the 
reference  ellipsoid  8h  and  the  time-correlated  barometric  altimeter  error  &ib. 

3.4.  Radar  Altimeter  Measurement  Model 

Accurate  and  precise  knowledge  of  altitude  is  vital  to  flying  and  landing  an  aircraft 
safely.  Thus,  an  additional  instrument,  a  radar  altimeter,  is  commonly  included  on 
aircraft  to  determine  altitude  further.  The  radar  altimeter  readings  hRad  are  modeled  as 
merely  the  sum  of  the  true  altitude  hj  and  zero-mean  WGN  in  the  measurement  Vnad  (note 
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the  sign  convention).  The  measurement  error  is  calculated  in  the  same  manner  as  with 
the  barometric  altimeter: 

^Rad  ~  ~  ^Rad 

=  [h,+Sh]-[h,-v,^,]  (67) 

=  Sh  +  v,„, 

Unlike  the  INS  and  barometric  altimeter  measurements,  there  are  no  time-correlated 
components  associated  with  the  radar  altimeter  measurements.  However,  the  radar 
altimeter  is  sensitive  to  the  relative  altitude  of  the  aircraft.  Radar  altimeters  give  better 
readings  when  closer  to  the  surface  terrain,  “above  ground  level”  (AGL),  than  at  higher 
altitudes.  Thus,  radar  altimeters  are  used  regularly  only  when  the  relative  altitude  drops 
below  3,000  ft  AGL.  The  precision  dependence  of  the  radar  altimeter  is  reflected  in  the 
covariance  of  the  measurement  noise  Rnad  as  a  function  of  the  true  AGL  altitude  /iagl* 
This  value  is  the  same  in  both  the  truth  and  filter  models  and  is  given  by  [13]: 

R^,=(0.0fhl^,+025)fe  (68) 

The  quantity  RRad  is  a  diagonal  element  of  the  measurement  noise  matrix  R. 

3.5.  Global  Positioning  System  Models 

GPS  is  made  up  of  a  constellation  of  satellites  sitting  in  six  orbital  planes  revolving 
about  the  earth  in  an  approximately  12-hour  orbit.  The  constellation  consists  of  a 
minimum  of  24  satellites  with  several  on-orbit  spares.  The  satellites  of  GPS  are  designed 
to  serve  as  a  radionavigation  system  to  provide  three-dimensional  positioning  information 
with  respect  to  the  earth.  Each  GPS  satellite  or  space  vehicle  (SV)  broadcasts 
information  about  its  own  position  in  orbit  at  a  specific  time.  GPS  receivers  on  earth 
search  for,  acquire,  and  track  signals  broadcast  from  four  or  more  SVs  and  process  them 
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to  determine  their  position  relative  to  the  earth.  Since  the  position  is  calculated  in  three 
dimensions,  at  least  four  signals  are  needed  to  produce  an  unambiguous  solution  for 
position.  The  signal  from  each  SV  contains  information  identifying  which  satellite  it  is, 
what  its  position  is  according  to  its  own  ephemeris,  what  time  this  signal  was  broadcast, 
and  other  data  specific  to  the  SV.  The  key  to  using  GPS  signals  for  positioning  is  to 
determine  the  difference  in  time  between  when  the  SV  sent  the  signal  and  when  the 
receiver  picked  up  the  signal.  Based  on  the  predetermined  velocity  of  radio  frequency 
(RF)  transmissions,  this  difference  in  time  can  be  translated  into  a  distance  between  the 
receiver  and  the  satellite.  This  quantity  is  known  as  the  pseudorange.  A  pseudorange 
differs  from  the  true  range  between  the  receiver  and  SV  in  that  it  has  several  errors 
embedded  into  it  due  to  timing  inaccuracies  and  RF  propagation  delays.  These  errors  are 
taken  into  account  when  determining  a  positioning  solution  for  the  receiver. 

3.5.1.  GPS  Truth  Models 

The  pseudoranges  calculated  from  the  receiver  to  each  GPS  satellite  are  subject  to 
three  primary  categories  of  errors:  errors  due  the  receiver’s  internal  clock,  errors  due  to 
atmospheric  effects,  and  errors  due  to  each  SV’s  equipment.  Each  error  has  the  effect  of 
increasing  or  decreasing  the  apparent  pseudorange  between  the  receiver  and  the 
transmitting  satellite.  Therefore,  these  errors  have  been  studied  widely  and  models  have 
been  constructed  to  describe  and  simulate  their  behavior  accurately  [25]. 

3.5.1.1.  GPS  Dynamics  Model 

The  GPS  truth  model  developed  has  been  used  with  success  by  several  research 
efforts  at  AFTT  [7, 13, 38, 56, 57].  The  30  states  used  by  the  model  are  listed  in 
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Appendix  B.  The  first  errors  considered  are  the  most  significant,  which  involve  the 
internal  clock  of  the  receiver.  Although  each  SV  in  the  GPS  constellation  possesses  a 
high-accuracy  atomic  clock  on-board,  most  receivers  use  less  expensive  clocks  that  do 
not  keep  time  to  the  same  degree.  The  difference  between  the  time  the  receiver  clock 
indicates  and  the  actual  time  is  represented  by  user’s  receiver  clock  bias  5Rucik  and  clock 
drift  rate  SDudk-  The  system  dynamics  equation  relating  these  is: 


'0 

r 

l^^Uclk 

A,.- 

0 

0 

(69) 


The  initial  values  of  the  state  estimates  and  corresponding  values  of  covariance  were 
chosen  to  be  consistent  with  past  studies  at  AFTT  [7, 13, 38,  56]: 
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These  errors  are  the  same  with  respect  to  each  satellite  since  they  are  produced  by  the 
receiver  equipment.  The  next  errors  considered  are  associated  with  a  particular  satellite, 
denoted  by  the  subscript  je  { 1, 2, 3, 4}.  Another  error  associated  with  GPS  receivers  is 
the  multi-path  error,  SRmpj,  which  is  caused  by  stray  signals  reflected  by  surfaces  near  the 
receiver’s  location. 

In  addition,  GPS  signals,  being  RF,  are  subject  to  atmospheric  interference  delaying 
the  transmission.  There  are  two  atmospheric  errors  modeled:  the  tropospheric  delay 
SRiropj  and  the  ionospheric  delay  5Rio„j.  The  three  errors  cited  above  are  all  treated  as 
first-order  Markov  processes  with  zero-mean  WGN  as  shown  by: 
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(72) 
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with  initial  state  covariances  given  by: 
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and  dynamics  noise  characteristics  given  by: 

£{Wg/.5(0}  =  0 


^{Wgp5  (OWgps  (t  +  T)}  =  QopsSir)  = 
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(75) 


Although  each  GPS  SV  has  a  high-accuracy  clock  on-board,  its  time  may  not  be  set 
exactly.  This  introduces  another  error  term  known  as  the  satellite  clock  bias  dRsdkj^  The 
last  errors  taken  into  account  are  the  line-of-sight  satellite  position  errors,  Sxsvj,  Sysvj-,  and 
Szsvj-  These  four  errors  are  modeled  thus: 


^Psclkj 

^^Sclkj 

SXgyj 

“  ®4x4 

Sx^YI 

^SVj 

^SVJ 

_  ^SVj  _ 

_  ^^SVj  _ 

with  initial  state  covariances  given  by: 


(76) 


Psc,uLos(to)  =  25feK,  (77) 

The  full  30-state  system  dynamics  equation  is  formed  by  taking  the  user’s  receiver 
clock  equations  and  augmenting  them  with  four  versions  of  Equations  (72)  and  (76),  for 
each  f'  satellite.  Appendix  B  provides  a  description  for  the  30  states. 
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3.5.I.2.  GPS  Measurement  Model 


The  principal  measurements  for  a  GPS  receiver  are  the  pseudoranges  popsj  to  each 
satellite.  This  can  be  expressed  as  the  sum  of  the  true  range  prj  plus  all  the  errors 
described  above  and  measurement  noise  Vprj: 

Pops  =  Prj  +*«  +''™; 

The  true  range  is  always  an  unknown  quantity,  thus  a  measurement  difference  is 
calculated  to  eliminate  it.  This  is  carried  out  by  computing  the  position  of  the  craft 
(user’s  receiver)  as  indicated  by  the  INS  pu  and  the  satellite’s  position  according  to  its 
ephemeris  ps.  The  difference,  INS  range,  is  given  by: 


Xjj 

ECEF 

Xs 

ECEF 

PlNS'^\Pu  P5I” 
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where  ECEF  indicates  that  all  elements  are  expressed  in  the  Earth-Centered,  Earth-Fixed 
(ECEF)  coordinate  reference  frame.  The  quantity  pm  also  can  be  expressed  as: 

Pm = V(%  -Xsf  +  (yu  -ysf + 

This  equation  for  pi^s  is  a  non-linear  function  of  the  craft  and  satellite  positions.  To 
create  a  linear  expression  from  this,  a  Taylor  series  expansion  is  calculated  about  a 
nominal  value  of  both  positions  with  the  terms  greater  than  first-order  neglected. 


This  yields  the  equation: 
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Substituting  Equation  (80)  into  Equation  (81)  and  evaluating  the  partial  derivatives 
produces: 


Pm  -  Pt  ( ( 


'  ys-yij 


(82) 


The  measurement  difference  equation  is  formed  by  taking  the  difference  between 


Equation  (82)  and  Equation  (78)  for  each  satellite: 

^  P/NSj  ~  Popsj 


teK-liSfK -(»)*» 


(83) 


~^PMPj  ~^^ropj  ~^PionJ  ~^PsclkJ  ~  ^^Uclk  +'^PR/ 

Thus,  the  values  of  the  receiver  position  errors,  5xu,  dyu,  and  dzu,  can  be  derived  from  the 
first  three  states,  Sdx,  59y,  and  36^,  of  the  truth  or  filter  model  through  an  orthogonal 
transformation  [7,  13]. 


3.5.2.  Differential  GPS  Truth  Models 

As  seen  in  the  previous  section,  GPS  pseudorange  data  contain  errors  from  many 
different  sources.  A  commonly  used  technique  to  eliminate  some  of  these  errors  is 
known  as  Differential  GPS  (DGPS).  To  apply  this  technique,  two  GPS  receivers  are 
required.  The  first  receiver,  known  as  the  reference  receiver,  is  located  at  a  site  whose 
position  is  well  surveyed  and  known  to  a  high  degree  of  accuracy.  The  reference  receiver 
is  usually  equipped  with  a  high-accuracy  clock  to  diminish  the  receiver  clock  bias  and 
drift.  By  comparing  its  GPS-computed  position  with  its  surveyed  position,  the  reference 
receiver  can  obtain  a  precise  estimate  of  the  errors  in  the  pseudoranges  to  each  visible 
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satellite.  For  example,  certain  errors,  such  as  the  error  in  the  position  of  the  SV  and 
satellite  clock  bias,  are  common  between  receivers  that  are  relatively  close  by  (within 
approximately  100  nmi).  These  errors  then  may  be  calibrated  and  transmitted  to  another 
remote  receiver,  e.g.,  on-board  an  aircraft,  in  the  vicinity  of  the  reference  receiver.  The 
transmitted  error  data  provide  differential  corrections  for  the  use  by  the  remote  receiver 
to  remove  them  from  its  pseudorange  calculation.  Using  DGPS,  the  truth  model  can  be 
reduced  to  26  states  by  disregarding  the  satellite  clock  biases  dRsdkj  for  each  satellite.  In 
addition,  there  is  difficulty  in  consistently  characterizing  the  multi-path  error,  SRmpj, 
since  it  is  highly  dependent  upon  the  location  of  the  receiver.  Therefore,  for  the 
simulations  performed  in  this  research,  SRmpj  is  lumped  into  the  measurement  noise,  V. 
Thus,  the  final  DGPS  truth  model  used  consists  of  22  states. 


3.5.2.I.  DGPS  Dynamics  Model 

The  system  dynamics  equation  for  the  receiver  clock  bias  and  drift  remains 
unchanged  as  in  Equations  (69),  (70),  and  (71).  Grouping  the  remaining  atmospheric  and 
satellite  position  errors,  the  system  dynamics  equation  becomes: 
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with  initial  state  covariances  given  by: 
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and  dynamics  noise  characteristics  given  by: 
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3.5.2.2.  DGPS  Measurement  Model 

The  DGPS  pseudorange  measurement  equation  is  identical  in  form  to  Equation  (78) 
less  the  range  errors  due  to  multi-path  and  satellite  clock  bias: 

PdGPSJ  ~  Ptj  ^^tropj  ^ionj  ■*"  ^Uclk  ~  ^ PRj  (^^) 

The  DGPS  measurement  difference  equation  is  formed  in  the  same  manner  as  previously 
seen.  In  this  case,  the  difference  between  Equation  (82)  and  Equation  (88)  is  taken  for 
each  satellite: 


^PRj  ~  PiNSJ  PdGPSJ 


3.5.3.  DGPS  Filter  Models 


The  simulations  conducted  in  this  work  utilize  DGPS  to  take  advantage  of  the  higher 
accuracy.  Of  all  the  errors  discussed,  the  predominant  errors  are  the  receiver  clock  bias 
SRucik  and  clock  drift  dRudk-  The  other  errors  are  combined  and  modeled  as  zero-mean 
WGN  resulting  in  this  system  dynamics  equation: 
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with  an  initial  covariance  given  by: 


Pt/dt(^o) 


9.0x10’^  ft^  0 

0  9.0xl0‘° 

sec^ , 


(91) 


To  get  the  filter  model  to  emulate  the  behavior  of  the  truth  model,  the  elements  of  the 
dynamics  noise  matrix  Q(t)  are  tuned  in  the  filter.  These  tuning  values  are  listed  in 
Appendix  C. 

The  measurement  equation  for  the  filter  is  similar  to  the  other  filter  equations  except 
several  states  are  omitted:  the  multi-path  effects,  atmospheric  effects,  satellite  clock 
biases,  and  satellite  position  errors  from  the  ephemerides.  With  only  the  error  in  the 
craft’s  (receiver’s)  position  and  receiver  clock  bias  and  drift  remaining,  the  filter 
measurement  equation  works  out  to  be; 

^PRj  ~  PlNSj  ~  PdGPSJ 

(92 

=  "(13  “  (ifej  "  ( w"  ^ 

The  strength  of  the  measurement  noise  is  adjusted  to  compensate  for  the  eliminated 
states.  The  four  diagonal  elements  of  the  measurement  noise  matrix  R  pertaining  to  the 
pseudorange  measurements  RprjJ  e  { 1,  2,  3, 4}  are  tuned  for  optimal  performance. 


54 


3.5.4.  GPS  Pseudolite  Models 


The  pseudorange  measurements  for  the  ground-based  pseudolite  are  treated  as  if  they 
were  received  from  a  fifth  GPS  satellite  (j  =  5),  with  three  exceptions  pertaining  to  the 
truth  model: 

1.  The  atmospheric  effects  on  pseudolite  transmissions  are  caused  only  by 
the  troposphere,  since  the  pseudolite  is  on  the  ground.  Thus,  there  is  no 
ionospheric  delay  error  term  dRions  present,  only  a  tropospheric  term 

^J^tropS- 

2.  The  pseudolite  is  assumed  to  be  located  at  a  surveyed  site,  eliminating  the 

.  uncertainty  in  the  satellite  position  from  its  ephemeris:  dxsvs,  Sysvs,  Szsvs,- 

3.  The  pseudolite  is  equipped  with  a  high-accuracy  clock,  thus  clock  bias 
SRscm  may  be  ignored. 

Otherwise,  for  the  filter  model,  there  are  no  differences  between  a  satellite  measurement 
and  a  pseudolite  measurement  that  need  to  be  taken  into  account. 


3.6.  Integrated  System  Models 

The  models  implemented  into  the  simulations  are  a  combination  of  those  described 
above.  The  state  vector  of  the  entire  system  is  constructed  so  the  first  13  states  of  the 
truth  state  vector  form  the  filter  state  vector.  The  remaining  states  for  the  truth  model  are 


augmented  afterwards,  in  the  following  format: 
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DGPS  Filter 
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1 1  INS  Filter  States 
2  DGPS  Filter  States 
Other  28  INS  States 
Other  21  DGPS  States 


(93) 
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The  elements  of  the  system  dynamics  matrix  F  and  the  dynamics  noise  strength  matrix  Q 
are  arranged  to  maintain  the  state  relationship  so  the  system  dynamics  equation  for  the 
truth  model  in  terms  of  the  filter  model  is: 

^^Fiher  ^Filler  ^Olher(ly2)  ^^Filter  ^  ^Filler 

^^Other_  ^Other(2,l)  ^Olher(2,2)  _  Other  _  1.^ Other 

and  the  noise  strength  is: 

£{Wr™,.(0wL*(^  +  'Z^)}  =  Qr™,/,(0<5(T)=  ^  8{t)  (95) 

L  ^  Mother  ^ 

All  of  the  dynamics  noise  for  each  state  is  considered  to  be  independent  of  (uncorrelated 
with)  each  other.  Therefore,  all  of  the  matrices  representing  the  strength  of  the  dynamics 
noise,  Qiruth,  QpiUer,  and  Qother,  are  diagonal. 

The  measurement  vector,  Sz  =  [  SzBar  Szpri  •  •  •  8zpr4  SzRad  Szprs  1^,  remains  the 
same  for  both  the  truth  model  and  the  filter  since  there  is  no  reduction  in  the  number  of 
measurement  sources  from  the  truth  model  to  the  filter.  Conversely,  the  measurement 
noise  covariance  matrix,  R,  does  have  different  values  between  the  truth  model  and  the 
filter  to  reflect  the  reduction  of  states  from  the  truth  model  to  the  filter.  In  addition,  the 
measurement  matrix,  H,  is  structured  to  fit  the  state  vector  convention.  Thus,  the 
measurement  equation  for  the  truth  model  in  terms  of  the  filter  is: 
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and  a  measurement  noise  strength  of: 
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where  Sy  is  the  Kronecker  delta  function. 


3.7.  Interference  Models 

The  environment  in  which  an  aircraft  travels  is  filled  with  many  sources  of  RF 
signals,  especially  near  the  ground.  Although  GPS  is  a  highly  accurate  system  for 
determining  positions  with  respect  to  the  earth,  it  processes  RF  signals  with  a  very  low 
carrier-to-noise  ratio,  -163  dBw  [25].  These  signals  may  be  subjected  to  interference 
from  various  sources  ranging  from  unintentional  spurious  emissions  to  malicious 
disruption  on  or  near  the  GPS  frequency  band.  In  these  simulations,  this  is  reflected  by 
an  increase  in  the  elements  of  the  measurement  noise  matrix  R(r,)  related  to  the 
pseudorange  measurements.  This  variation  in  R(t,)  represents  a  change  in  the  parameter 
structure  of  the  system  a(t,).  Accounting  for  this  change  while  generating  the  appropriate 
state  estimates  is  the  role  of  the  M^AE.  The  variable  elements  RpRjiti)  of  R(t/)  will  have  a 
finite  range,  i.e.,  parameter  space,  since  beyond  a  certain  level  of  interference,  GPS  signal 
are  effectively  unusable.  The  MMAE  is  designed  to  have  filters  spanning  this  finite 
parameter  space,  and  to  be  optimized  for  estimating  the  actual  value  of  R(t,)  based  on  the 
measurements  available  up  to  the  time  in  question,  Z(t,)  =  [z  (t,)  Z  (t/.;)  •  •  •  z  (to)]  .  No 
changes  in  the  measurement  noise  for  the  altimeters,  Rsariti)  and  RRadiU),  are  taken  in 
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account.  The  barometric  altimeter  is  affected  only  by  the  atmospheric  pressure  and  the 
radar  altimeter  has  a  narrow  field  of  reception  directed  downward.  In  addition,  the  INS  is 
a  self-contained  piece  of  equipment  influenced  only  by  the  linear  and  angular 
accelerations  of  the  craft.  With  a  feed-forward  configuration  in  place,  the  only 
measurement  the  INS  receives  is  from  the  barometric  altimeter  incorporated  within  itself. 

3.8.  Simulation  Software 

The  FORTRAN-based  software.  Multimode  Simulation  for  Optimal  Filter  Evaluation 
(MSOFE),  is  a  program  constructed  for  conducting  simulations  of  systems  that  employ 
Kalman  filtering  [36].  The  structure  of  MSOFE  allows  it  to  be  directly  applied  to  a 
variety  of  problems  requiring  optimal  state  estimation  with  a  minimal  amount  of  software 
refitting.  The  program  has  two  principal  simulation  modes  for  assisting  in  evaluating  the 
performance  of  the  system: 

1.  Monte  Carlo  simulation  generates  multiple-sample  time  histories  of  a 
system’s  truth  states,  filter  states,  and  filter  estimation  errors,  including 
non-linear  effects,  using  linear  or  extended  Kalman  filtering. 

2.  Covariance  simulation  generates  time  histories  of  the  second-order 
statistics,  i.e.,  values  of  covariance,  of  a  system’s  truth  states,  filter 
states,  and  filter  estimation  errors  using  only  linear  or  linearized 
models. 

The  two  modes  are  complementary.  Covariance  simulation  can  generate  the 
performance  statistics  of  a  filter  in  a  single  pass.  Monte  Carlo  simulation  requires 
multiple  sample  runs  before  meaningful  statistics  can  be  computed  for  a  scenario. 
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Conversely,  Monte  Carlo  simulation  can  handle  non-linear  processes,  which  the 
covariance  mode  cannot. 

An  extension  of  MSOFE,  Multiple  Model  Simulation  for  Optimal  Filter  Evaluation 
(MMSOFE),  developed  by  Nielsen  [39],  supports  the  analysis  of  systems  using  multiple 
model  adaptive  estimation.  MMSOFE  is  a  modified  version  of  MSOFE  using  similar 
core  code  with  modifications  reflecting  the  use  of  multiple  Kalman  filters.  It  is  designed 
to  conduct  the  propagation  and  update  of  several  filters  simultaneously.  While 
performing  the  filter  operations,  it  performs  the  calculations  for  hypothesis  conditional 
probability  and  appropriate  blending  for  MMAE  and  other  multiple  model  algorithms.  In 
this  work,  MMSOFE  is  used  in  the  Monte  Carlo  mode  for  all  simulations. 

3.9.  Parameter  Estimation 

The  navigation  system  models  are  implemented  using  the  M^AE  architecture  with  an 
MMAE  consisting  of  five  elemental  filters,  delivering  parameter  updates  to  a  state 
estimator.  Five  was  chosen  to  give  a  finer  discretization  of  the  parameter  space  than  in 
the  previous  research  by  White  [56],  but  to  keep  the  simulation  time  needed  to  a 
manageable  level.  The  estimated  parameter  is  the  GPS  pseudorange  measurement  noise 
strength,  RprjJ  e  { 1, 2, 3, 4, 5}.  This  noise  is  assumed  to  affect  all  pseudorange 
measurements  uniformly,  thus  Rppi  =  Rppj  V  i  ^j.  Although  the  values  of  Rpr  (j  subscript 
dropped  for  simplicity)  are  meant  to  portray  an  increase  in  the  strength  of  the  RF 
interference  noise  on  the  GPS  LI  frequency,  it  is  not  an  accurate  representation  of  an 
increase  in  the  noise  power.  Still,  varying  Rpr  serves  to  demonstrate  a  meaningful 
change  in  a  parameter. 
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The  greatest  admissible  value  of  pseudorange  measurement  noise  variance  is  set  as 

2,000  times  the  nominal  value  of  Rpr  without  RF  interference  (Rpro)  [56].  Any  value 

greater  than  this  is  assumed  to  render  GPS  signals  useless.  Thus,  the  finite  parameter 

space,  K,  can  be  defined  to  encompass  the  real  numbers  from  1  to  2,000  and  the 

estimated  parameter  vector,  a,  is  actually  a  scalar  quantity,  a  €  S .  The  value  of  a  can  be 

treated  then  as  a  multiplier  to  Rpro,  which  in  past  research  at  AFTT  was  taken  to  be  9  (ft ) 

[38].  With  K  and  a  defined,  the  placement  of  the  elemental  MMAE  filters  can  be 

determined.  Using  the  algorithm  developed  by  Sheldon  [48]  described  in  Section  2.4.1., 

the  values  of  a  (instead  of  a,  since  this  quantity  is  not  stochastic)  for  each  of  the  five 

filters  was  calculated  as  follows: 

a,  =  1.0 
02  =  214.4 
O3  =  610.7 
O4  =  1,064.3 
Oj  =1,629.4 

In  addition,  a  constrained  optimization  was  implemented  into  the  enhanced  Sheldon 
algorithm  to  fix  one  of  the  filters  to  the  nominal  parameter  level,  a\. 

During  each  Monte  Carlo  run,  three  programs  are  used.  First,  the  program  MMSOFE 
runs  the  MMAE  filters  to  calculate  its  best  estimate  of  the  pseudorange  measurement 
noise  strength  factor,  d(ti)  =  £{a(t/)  |  Z(r/)  =  Z,  },  at  each  time  interval  over  the  entire 
flight  trajectory  in  question.  This  value  of  d  is  then  appended  to  the  flight  trajectory  data 
using  the  program  APPEND  and  passed  on  to  the  state  estimator  to  be  incorporated  into 
its  measurement  noise  matrix,  R,  by  letting  Rrrj  =  dRppo.  Next,  the  state  estimator 
calculates  its  best  estimate  of  the  states,  x  ,  from  which  the  best  estimate  of  the  error  in 
the  aircraft’s  position  is  calculated  at  each  time  interval.  Like  the  MMAE,  the  state 
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estimator  is  implemented  into  MMSOFE  (a  variation  called  SMSOFE)  with  certain 
variations,  e.g.,  the  number  of  filters  is  set  to  one.  SMSOFE  has  the  same  basic  structure 
as  MMSOFE,  except  it  has  only  one  filter  that  accepts  a  parameter  update  generated  by 
the  MMAE  (implemented  in  MMSOFE).  After  SMSOFE  is  run,  the  error  in  the  aircraft’s 
position  (east,  north,  and  vertical  directions)  over  the  trajectory  is  recorded.  Ten  Monte 
Carlo  runs  are  performed  to  obtain  sufficient  data  for  a  meaningful  analysis.  Finally, 
these  data  are  processed  using  programs  written  in  Matlab®  [52]  to  create  plots  of  the 
errors  over  the  time  of  the  trajectory. 

3.10.  Summary 

The  models  described  in  this  chapter  are  employed  for  simulating  an  integrated 
GPS/INS  aircraft  navigation  system  augmented  with  a  radar  altimeter  and  a  signal 
available  from  a  GPS  pseudolite.  A  reduced-order  filter  model  is  constructed  using  a 
subset  of  the  states  from  the  system  truth  model.  This  is  done  to  model  the  difference 
between  the  predictions  an  on-board  navigation  computer  (filter  model)  would  make  of 
the  position  of  the  aircraft  against  actual  behavior  of  the  craft  (tmth  model).  In  Chapter 
4,  the  performance  of  the  system  is  analyzed  through  various  levels  of  interference  to 
determine  its  suitability  for  a  PLS. 
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4.  Analysis  of  Results 


This  chapter  covers  the  simulations  conducted  and  provides  interpretations  of  the 
results.  Eight  test  cases  are  presented,  depicting  two  navigation  system  configurations 
and  four  GPS  interference  scenarios.  The  M^AE  algorithm  is  implemented  for  each  case, 
where  the  MMAE  is  set  up  to  determine  the  variable  parameter  in  the  presence  of 
interference  and  the  state  estimator  provides  the  best  estimate  of  the  position  of  the 
aircraft. 

4.1.  Flight  Profile 

The  simulations  were  performed  using  the  MMSOFE  program  working  with  a  flight 
trajectory  created  by  PROFGEN  [35].  The  flight  profile  created  represents  a  KC-135 
(Boeing  707)  aircraft  flying  for  3,912.5  seconds  from  take-off  to  landing  touch-down  at 
runway  23R  (heading  232°)  at  Patterson  Field,  Wright-Patterson  AFB.  The  portion  of  the 
flight  profile  examined  is  the  period  just  before  landing,  from  3,810  to  3,910  seconds, 
when  the  aircraft  is  on  its  3°  descent  glideslope  with  a  steady  velocity  of  132  knots  [13]. 
The  local  terrain  is  assumed  to  be  “flat”,  at  an  altitude  of  825  ft  from  mean  sea  level 
(MSL).  The  coordinates  of  the  flight  path  are  shown  in  Figure  10. 
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Actual  Aircraft  Trajectory 


Time  (sec) 

Figure  10.  Actual  Trajectory  of  the  Simulated  Landing  Aircraft 


4.2.  GPS  Data 

To  determine  the  pseudorange  measurements,  actual  GPS  almanac  data  were 
incorporated  into  the  simulations.  The  almanac  data  were  obtained  from  the  Coast  Guard 
Bulletin  Board  Services  (CGBBS).  A  program  called  SEM3.6  [13]  processed  this 
information  to  select  the  four  best  GPS  satellites  of  those  visible,  based  on  the  best 
geometry,  i.e.,  geometric  dilution-of-position  (GDOP).  An  arbitrary  take-off  time  and 
date,  04:00  UTC  (08:00  EDT),  21  May  1994,  and  a  mask  angle  of  5°  above  the  horizon, 
were  used  to  determine  which  GPS  satellites  would  be  visible  for  the  flight  trajectory. 
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4.3.  Navigation  System  Configurations 

The  navigation  systems  modeled  in  this  study  follows  the  work  of  previous  [7, 13, 38, 
56]  and  current  [53,  57]  research  at  AFTT.  The  core  of  the  system  is  an  INS  with  a 
barometric  altimeter  built-in  for  vertical  channel  aiding.  To  compare  the  effects  of  the 
quality  of  the  INS  used,  two  models  were  used,  an  INS  with  a  drift  rate  of  0.4  nmi/hr 
(CEP)  vs.  a  lower  quality  INS  with  a  drift  rate  of  4.0  nmi/hr  (CEP).  The  positioning 
ability  of  the  INS  is  enhanced  with  pseudorange  measurements  provided  by  a  GPS 
receiver.  The  receiver  accepts  signals  from  the  four  best  GPS  satellites  visible  and  from  a 
ground-based  GPS  pseudolite.  The  pseudolite  is  included  to  improve  the  vertical 
positioning  of  GPS.  This  is  reflected  in  a  lower  (better)  value  of  GDOP;  2.63  without  the 
pseudolite  vs.  2.50  using  the  pseudolite  [25].  Finally,  a  radar  altimeter  is  incorporated  to 
improve  altitude  readings  further;  this  is  a  navigation  aid  found  in  most  military  and 
commercial  aircraft  for  use  in  landing. 

4.4.  GPS  RF  Interference  Scenarios 

The  navigation  systems  modeled  in  this  study  follows  the  work  of  previous  research 
by  Gray,  Britton,  and  White.  Gray  examined  different  combinations  of  navigation  aids, 
first  introducing  different  qualities  of  INS:  0.4  nmi/hr,  2.0  nmi/hr,  and  4.0  nmi/hr  drift 
rates.  Britton  followed  this  by  incorporating  DGPS  in  place  of  ordinary  GPS 
pseudorange  measurements.  White  applied  RF  interference  and  GPS  spoofing  to  some  of 
the  past  test  cases  to  assess  the  performance  of  an  MMAE  to  provide  accurate  state 
estimates. 
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The  simulations  performed  for  this  effort  use  four  different  schemes  of  RF 
interference  to  examine  the  performance  of  the  M^AE  design  described  in  previous 
chapters.  The  first  has  no  interference  to  establish  a  nominal  reference.  The  second  uses 
the  same  scheme  applied  by  White  in  his  research,  which  consists  of  varying  the 
parameter,  a,  at  different  levels  spaced  as  shown  in  the  top  plot  of  Figure  11.  The  third 
scheme  attempts  to  portray  a  situation  where  the  interference  increases  as  the  aircraft 
approaches  the  touch-down  point.  Here,  the  source  is  assumed  to  be  in  the  flight  path  of 
the  aircraft.  The  value  of  a  is  increased  in  inverse  proportion  to  the  square  of  the  distance 
to  the  interference  source.  In  the  fourth  scheme,  the  value  of  a  is  decreased  in  the  same 
manner  to  depict  flying  away  from  an  RF  interference  source.  The  four  RF  interference 
scenarios  are  combined  with  the  two  navigation  configurations  to  form  eight  test  cases. 
Table  3  summarizes  the  test  cases  and  Figure  1 1  shows  the  RF  interference  schemes 
where  the  parameter  a  varies. 


Table  3.  Simulation  Test  Cases 


Case 

Drift  Rating  of 
INS  Used 

RF  Interference 
Scheme 

I 

0.4  nmi/hr 

Nominal 

n 

4.0  nmi/hr 

Nominal 

m 

0.4  nmi/hr 

White’s  (“Exponential”) 

IV 

4.0  nmi/hr 

White’s  (“Exponential”) 

V 

0.4  nmi/hr 

Increasing 

VI 

4.0  nmi/hr 

Increasing 

VII 

0.4  nmi/hr 

Decreasing 

vm 

4.0  nmi/hr 

Decreasing 
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Exponential  Noise  Levels  Scenario:  Cases  111  &  IV 
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Decreasing  Noise  Level  Scenario:  Cases  VI  &  VIII 


Increasing  Noise  Level  Scenario:  Cases  V  &  VI 


Figure  11.  GPS  Interference  Scenarios 

4.5.  PLS  Comparison  Criteria 

When  the  aircraft  is  on  final  approach,  the  FAA ILS  criteria  [12]  are  applied  at 
specific  decision  heights  (DH).  When  the  aircraft  reaches  a  particular  DH,  the 
uncertainty  in  the  horizontal  and  vertical  must  be  no  greater  than  the  values  specified  by 
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that  category.  These  values  and  the  time  in  the  trajectory  when  the  aircraft  in  these 
simulations  reaches  the  DH  points  are  summarized  in  Table  4. 


Table  4.  FAA ILS  Precision  Approach  Requirements  (±la) 


Category 

Decision  Height 
(AGL) 

Occurrence  in 
Flight  Trajectory 

Horizontal 

Accuracy 

Vertical 

Accuracy 

I 

200  ft 

3,894  sec 

±28.1  ft 

±  6.8  ft 

II 

100  ft 

3,902  sec 

±  8.6  ft 

±  2.8  ft 

III 

50  ft 

3,907  sec 

±  6.8  ft 

±1.0  ft 

Two  types  of  data  plots  are  displayed  in  the  next  sections  covering  each  test  case. 

The  first  type  includes  six  graphs.  The  top  graph  is  a  plot  of  the  MMAE’s  estimate  of  the 
parameter,  d  over  the  simulation  period.  The  solid  line  indicates  the  mean  value  of  d  over 
ten  Monte  Carlo  runs,  while  a  dashed  line  above  and  below  show  the  span  of  one 
standard  deviation  in  d  (±(Ta)  over  the  ten  Monte  Carlo  runs.  The  heavier  solid  flat  lines 
represent  the  true  values  of  a.  The  horizontal  dotted  lines  show  the  values  of  a  where 
each  elemental  filter  is  placed,  i.e.,  a*,  e  { 1,  2,  3, 4,  5}.  The  five  graphs  below  the 
estimated  parameter  plot  illustrate  how  the  probability  assigned  to  each  elemental  filter, 
Pk,  varies  over  time.  Again,  a  solid  line  represents  the  mean  value  of  the  quantity  in 
question,  pk,  and  the  dashed  lines  represent  the  span  of  one  standard  deviation  of  pk  (±<Tp) 
over  ten  Monte  Carlo  runs.  Figure  12  shows  a  legend  for  these  plots. 

The  second  type  of  figure  shows  three  graphs.  Each  one  displays  the  error  in  the 
position  of  the  aircraft  as  it  varies  during  the  landing  period.  The  top  graphs  shows  the 
variation  in  error  in  the  east  horizontal  (longitude)  direction,  the  middle  graph  shows  the 
same  for  the  north  horizontal  (latitude)  direction,  and  the  bottom  graphs  shows  this  for 
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the  vertical  (altitude)  direction.  The  solid  line  in  the  center  of  each  plot  represents  the 
mean  position  error  over  ten  Monte  Carlo  runs.  The  dotted  lines  above  and  below  it 
represent  the  span  of  one  standard  deviation  of  the  position  error  {±(7True)  over  the  ten 
Monte  Carlo  runs.  This  is  the  actual  level  of  estimation  performance  encountered  in  the 
simulations.  The  dashed  lines  above  and  below  these  represent  the  span  of  one  standard 
deviation  computed  by  the  filter  (±<yFiiter)-  According  to  the  filter  model,  this  is  the  level 
of  uncertainty  expected  under  the  conditions  of  the  simulation.  In  addition,  the  filter- 
computed  standard  deviation,  Ofuter,  is  shown  added  to  and  subtracted  from  a  value  of 
zero,  since  the  filter  assumed  the  errors  to  be  a  zero-mean  process.  Figure  13  shows  a 
legend  for  these  plots. 

Filter-Computed  Parameter/Probability 

Standard  Deviation  over  10  Monte  Carlo  runs 

Parameter  Value  of  an  Elemental  MMAE  Filter 

Figure  12.  Parameter  /  Estimated  Probability  Plot  Legend 

Position  Error  over  10  Monte  Carlo  runs 

Standard  Deviation  over  10  Monte  Carlo  runs 

Filter-Computed  Standard  Deviation 


Figure  13.  Position  Error  Plot  Legend 
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4.6.  Nominal  Interference  Scenario 

The  two  cases  that  follow  are  set  up  to  examine  the  performance  of  the  M^AE  model 
in  the  absence  of  RF  interference.  The  parameter  a  is  set  to  1  to  represent  a  condition 
where  the  pseudorange  measurement  noise  variance  is  equal  to  the  nominal  value,  i.e., 
Rprj  =  Rpro.  ;■  e  { 1, 2,  3, 4, 5 }  (see  Section  3.9).  Although  a  =  1,  the  MMAE  estimated 
parameter  delivered  to  the  state  estimator,  d,  is  a  higher  value.  This  is  due  to  the 
minimum  probabilities  assigned  to  each  elemental  filter,  Pmin*  The  MMAE  implemented 
in  these  simulations  have  pnun  =  0.001.  Thus,  assuming  the  a\  elemental  filter  receives  all 
of  the  probability  weight  less  the  minimum  assigned  to  the  other  filters,  the  estimated 
parameter  used  by  the  state  estimator  is  computed  by: 

a  =  =0.996a,+0.00lXa*  =4.514797  (98) 

k=l  k=2 

Although  this  value  is  over  four  times  greater  than  the  nominal,  it  is  the  best  the  MMAE 
can  do  given  the  pmin  specified.  Given  the  nature  of  this  problem,  it  recommended  that 
the  blending  could  be  limited  to  remove  the  effects  caused  by  pmin  (see  Section  5.3). 


4.6.1.  Case  I 

This  case  covers  the  best  scenario  of  the  ones  examined,  using  the  INS  with  the  lower 
drift  rate,  0.4  nmi/hr,  and  having  no  added  measurement  noise.  Figure  14  shows  how  the 
estimated  parameter,  d,  and  the  elementary  filter  probability  weights,  pk,  vary  over  time. 
Without  RF  interference  present,  the  MMAE  parameter  estimator  decides  very  quickly 
that  the  ai  elemental  filter  is  the  best  model  for  the  system.  Since  a  never  changes 
throughout  the  scenario,  d  remains  unchanged  once  the  MMAE  determines  the  a\ 
elemental  filter  has  the  best  representation  of  the  situation.  The  initial  spike  in  the  plot  of 


69 


Parameter  Value 


d  is  due  to  the  fact  that  each  elemental  filter  has  the  same  initial  probability  assigned  to  it, 
as  shown  in  Figure  14.  At  3,81 1  seconds,  the  az  elemental  filter  tries  to  grab  all  of  the 
probability  weight,  but  it  is  not  the  best  representation  of  the  current  environment.  The 
MMAE  adapts  quickly  once  the  ai  elemental  filter  residuals  look  the  best,  reflecting  the 
current  operating  environment. 

Figure  15  shows  the  output  of  the  state  estimator.  The  standard  deviation  of  the 
errors  in  all  three  cardinal  directions  remains  under  5  ft  for  nearly  the  entire  period  in 
question.  The  filter-computed  standard  deviation,  GFUter,  plots  (dashed  lines)  start  at  a 
conservatively  high  value.  It  converges  to  a  mean  steady-state  value  of  1.02  ft  in  the 
east(/west)  direction  and  1.30  ft  in  the  north  (/south)  direction  20  seconds  into  the 
trajectory.  The  value  of  ^Filter  in  the  vertical  direction  quickly  converges  to  a  value  of 
4.02  ft  within  5  seconds  and  decreases  steadily,  as  the  aircraft  approaches  ground  level. 
This  to  due  to  the  improved  measurement  estimates  in  the  vertical  direction  (i.e.,  altitude) 
due  to  the  radar  altimeter.  Equation  (68)  shows  that  the  variance  of  the  noise  in  the  radar 
altimeter  measurement,  Rnad,  is  a  function  of  the  altitude  AGL,  Hagl- 

(O.Of /lie. +0.25)  ft^ 

Thus,  the  closer  the  aircraft  is  to  the  ground,  the  better  the  radar  altimeter  readings 
become.  This  improvement  was  demonstrated  by  the  research  of  Gray  [13]  and  Britton 
[7],  where  the  value  of  CFUter  converged  to  a  steady-state  value  and  remained  unchanged 
when  radar  altimeter  measurements  were  not  incorporated.  Their  research  showed  the 
same  tjqje  of  decrease  in  OFUter  in  the  vertical  direction  seen  in  Figure  15  when  radar 
altimeter  measurements  were  included. 
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The  standard  deviation  of  position  errors  over  ten  Monte  Carlo  runs,  arrue,  (dotted 
lines)  is  consistently  less  than  aputer.  The  indicates  that  the  state  estimator  overestimates 
the  error  in  position,  a  conservative  tuning  philosophy.  The  mean  of  the  position  errors, 
even  over  ten  Monte  Carlo  runs  (solid  line),  does  not  remain  exactly  at  zero,  ideally,  it 
should  be  (this  was  seen  even  over  20  Monte  Carlos  runs;  ten  runs  were  chosen  for  each 
case  due  to  the  length  of  computer  simulation  time  required).  Still,  over  time,  the 
position  errors  do  stay  close  to  zero.  In  addition,  as  long  as  the  mean  position  errors  ± 
(jTrue  remain  close  to  the  (yputer  bounds,  the  state  estimator  filter  provides  good  estimates, 
as  shown  by  Figure  15. 

The  critical  tests  come  at  the  DHs  of  the  FAA ILS  categories.  Tables  5-7  show  the 
values  of  Cputer  and  (JTrue  at  each  DH  for  the  three  position  errors.  To  assess  horizontal 
accuracy,  the  values  of  OpiUer  in  the  east  and  north  directions  are  compared  to  the  FAA 
requirements.  In  Case  I,  the  values  of  both  (Jputer  and  Gprue  satisfy  the  Category  HI 
requirement,  ±  6.8  ft.  In  the  vertical  direction,  (Jputer  and  (JTrue  also  satisfy  the  Category 
TTT  requirement,  being  ±1.0  ft.  The  performance  of  Case  I  simulates  the  best  situation  of 
all  the  test  cases  considered.  Thus,  it  serves  as  a  reference  to  determine  how  a  lower 
quality  INS  or  RF  interference  degrades  the  performance  of  the  navigation  system. 

Table  5.  Case  I:  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

(Jtrue  (ft) 

I 

200 

1.034 

0.263 

n 

100 

1.034 

0.142 

m  1 

50 

1.034  ^ 

0.302 
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Table  6.  Case  I:  North  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

OnUer  (ft) 

Olrue  (ft) 

I 

200 

1.345 

0.342 

n 

100 

1.361 

0.233 

III  n 

50 

1.369 

0.530 

Table  7.  Case  I:  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

OFiUer  (ft) 

C^Tue  (ft) 

I 

200 

1.814 

1.214 

II 

100 

1.145 

0.711 

m 

50  1 

0.709 

0.798 

4.6.2.  Case  II 

This  case  is  identical  to  Case  I,  except  an  INS  with  a  drift  rate  of  4.0  nmi/hr  is 
implemented  in  the  navigation  system  model.  Figure  16  shows  essentially  identical 
results  to  Figure  14,  for  d  and  pk.  This  follows,  since  an  INS  is  not  affected  by  any 
external  measurements  and,  the  drift  rate  of  the  INS  is  not  a  function  of  the  measurement 

noise  (since  it  is  not  being  updated  by  the  filter). 

Conversely,  this  does  not  apply  to  the  output  of  the  state  estimator,  as  shown  in 
Figure  17.  While  the  position  errors  in  all  three  directions  remain  under  5  ft,  the  plots  of 
(Tpiiter  show  that  its  steady-state  values  are  greater  in  the  horizontal  directions;  1.83  ft, 
east,  and  1.99  ft,  north.  The  values  of  <7r™«  are  still  significantly  smaller  than  OFUteA 
however,  as  the  position  errors  wander  about  zero,  they  fit  closely  within  the  ^Filter 
bounds. 

Tables  8-10  show  how  OFUter  and  OTme  compare  with  the  ILS  categories.  In  both  the 
horizontal  and  vertical  directions,  the  values  of  OFiuer  and  Ojrue  satisfy  the  Category  III 
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Table  8.  Case  II:  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

OTrue  (ft) 

I 

200 

1.834 

0.593 

n  n 

100 

1.829 

0.348 

m 

50  ^ 

1.826 

0.457 

Table  9.  Case  II:  North  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

(^Filter  (ft) 

OTrue  (ft) 

I 

200 

1.999 

0.723 

n 

100 

1.999 

0.354 

m  n 

50 

1.997  n 

0.808 

Table  10.  Case  II:  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

(^Filter  (ft) 

OTrue  (ft) 

I 

200 

1.914 

1.836 

n 

100 

1.246  n 

0.900 

m  n 

50 

0.792 

0.824 

requirements.  This  shows  that  an  INS  with  a  greater  drift  rate  may  be  used  to  met  the 
requirements  with  nominal  GPS  pseudorange  measurement  noise  present.  In  addition, 
this  case  serves  as  a  reference  for  Cases  IV,  VI,  and  Vm  to  determine  how  increased 
level  of  pseudorange  measurement  noise  (i.e.,  increase  in  a,  hence  Rprj)  affects  the 
performance  of  this  navigation  system  (with  the  4.0  nmi/hr  INS). 


4.7.  White’s  (“Exponential”)  Noise  Levels  Scenarios 

The  next  two  cases  examine  the  system  performance  against  a  scenario  used  in 
White’s  research  [56].  In  this  scenario,  the  value  of  the  parameter  a  varies  at  levels  of  1, 
100, 200, 400,  800, 1600,  and  2000  during  the  100-second  interval  under  investigation;  as 
shown  in  the  top  plot  of  Figure  1 1 .  White  used  an  MMAE  with  three  elemental  filters 
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with  parameters  spaced  at:  a,  =  1,  ^2  =  200,  and  as  =  2,000.  The  position  errors  (states) 
calculated  by  each  filter  were  multiplied  by  its  corresponding  probability  weight,  pk,  k  e 
{ 1, 2, 3},  and  then  summed  together  to  form  a  blended  output. 

In  the  cases  performed  in  this  research,  blending  is  done  to  calculate  d  to  pass  on  to 
the  state  estimator,  which  determines  the  position  errors.  The  difference  between  the 
position  errors  of  a  blended  output  of  the  MMAE  and  those  of  the  state  estimator  are 
examined  in  this  case  to  provide  an  initial  baseline.  However,  it  is  important  to  note  that 
the  MMAE  is  tuned  for  parameter  estimation  in  the  M^AE  architecture. 


4.7.1.  Casein 

This  case  involves  the  navigation  system  using  the  0.4  nmi/hr  drift  rate  INS  with  the 
parameter  a  (thus,  Rprj)  at  levels  researched  by  White  [56].  Figure  18  shows  how  d  and 
Pk  vary  over  the  scenario.  In  the  top  plot,  the  actual  value  of  the  parameter  a  is  shown  by 
the  heavy  solid  lines.  The  MMAE  estimate,  d,  adapts  its  value  whenever  a  changes. 
During  the  first  10  seconds,  the  MMAE  estimates  d  =  4.51  as  in  the  nominal  scenario 
when  a  =  1.  When  a  becomes  100  in  the  next  10  seconds  (3,820  -  3,830  sec),  d  settles 
on  a  mean  value  of  216.9  within  2  seconds.  This  is  very  close  to  the  value  of  ai  =  214.4, 
indicating  the  MMAE  weighted  the  a2  elemental  filter  with  nearly  all  of  the  probability 
weight  (p2  =  0.996)  and  judged  it  to  be  the  most  accurate  model  for  the  current  operating 
environment.  Thus,  minimal  blending  occurred  during  this  period.  Also,  the  standard 
deviation  of  d,  a^,  remains  below  0.1,  showing  that  the  MMAE  was  very  consistent  m 
assigning  d  =  216.9.  This  can  be  seen  by  the  five  probability  plots;  p2  receives  almost  all 
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of  the  probability  (only  0.996,  though,  due  to  pmm  of  the  other  filters)  very  quickly  while 
the  elemental  filters  have  probabilities  near  zero. 

For  the  period  when  a  =  200  (3,830  -  3,840  sec),  the  mean  of  a  stays  at  216.9,  which 
is  appropriate  since  that  value  is  close  to  02.  When  a  changes  level  to  400  and  800  (3,840 
-  3,850  sec  and  3,850  -  3,860  sec,  respectively),  the  MMAE  starts  to  show  some 
blending.  The  value  of  d  does  not  quite  reach  a,  though  it  does  increase.  The  probability 
plots  show  the  probability  weighting  being  distributed  primarily  between  pi  and  ps-  The 
increase  in  the  standard  deviations  of  a  and  pk  indicate  the  MMAE  may  favor  one 
elemental  filter  over  the  others  in  a  single  Monte  Carlo  run.  Still,  taking  the  mean  over 
the  ten  Monte  Carlo  runs,  as  shown  in  the  plots  of  Figure  18,  it  averages  out  to  show  that 
probability  blending  taking  place. 

As  3  reaches  the  levels  of  1,600  and  2,000  (3,860  —  3,870  sec  and  3,870  —  3,880  sec, 
respectively),  the  MMAE  takes  all  of  the  period  when  a  =  1,600  to  push  d  near  that  value, 
and  d  hits  a  ceiling  at  1,612.5  when  a  =  2,000.  This  ceiling  is  due  to  the  highest  value  of 
Uk  for  an  elemental  filter,  which  is  as  =  1,629.4.  The  parameter  estimate,  d,  can  never  be 
greater  than  as,  and  its  upper  limit  is  slightly  less  due  to  the  pmin  that  must  be  assigned  to 
the  other  filters.  Even  though  the  MMAE  cannot  estimate  a  =  2,000  exactly  in  the 
configuration  used  in  this  research,  the  M^AE’s  state  estimate  performance  is  still 
effective  as  shown  later. 

For  the  remainder  of  the  scenario,  when  a  decreases  from  800  to  1,  the  mean  of  d 
follows,  but  with  a  lag  of  approximate  5  seconds.  This  is  a  typical  aspect  of  fixed-bank 
MMAEs.  The  delay  occurs  since  the  majority  of  the  probability  weight  resides  with  the 
filter  that  has  the  higher  ak  value  and  consequently  smaller  L^ft,)  =  rj value. 
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This  gives  the  MMAE  the  false  impression  that  it  has  the  best  parameter  estimation 
performance  [33].  However,  the  MMAE  quickly  responds  once  the  residuals  indicate 
this  performance  is  incorrect.  Subsequently  the  estimation  performance  improves  when 
the  elemental  filter  with  the  best  estimate  of  a  is  given  the  majority  of  the  probability 
weight. 

Figure  19  shows  the  blended  output  of  the  position  errors  from  the  MMAE  (tuned  for 
parameter  estimation).  Although  the  MMAE  used  in  these  cases  is  not  tuned  for  state 
estimation,  the  output  is  examined  to  give  a  link  to  the  past  research  of  White  [56],  which 
involved  an  MMAE  used  for  state  estimation.  These  plots  follow  the  results  of  White’s 
[56]  research,  using  his  navigational  set-up  and  scenario  (cf.  Nav  Case  1  in  [56]).  The 
minpr  differences  between  the  plots  here  and  in  [56]  can  be  attributed  to  using  a  different 
set-up  of  the  elemental  filters  in  the  MMAE.  Blending  does  provide  considerably  better 
estimates  of  the  position  errors  than  a  single  non-adaptive  Kalman  filter  could.  Figure  20 
shows  the  position  error  plots  from  the  output  of  the  ai  elemental  filter  (note  the 
difference  in  scale).  This  filter  measurement  noise  parameter  is  tuned  for  performance 
under  nominal  conditions.  When  the  measurement  noise  changes,  as  in  this  case,  its 
performance  drastically  suffers.  Improving  performance  under  changing  parameters  is 
the  thrust  for  developing  adaptive  filtering  techniques  such  the  M  AE. 

Figure  21  shows  the  positions  errors  calculated  by  the  state  estimator  within  the 
M^AE.  The  horizontal  (east  and  north)  position  error  estimates  improve  when  using  a 
filter  provided  with  an  estimate  of  the  uncertain  parameter  (the  state  estimator)  over  using 
ordinary  MMAE  blending  techniques.  As  mentioned  before,  the  MMAE  used  in  these 
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Figure  21.  Case  III:  State  Estimator  Position  Errors 
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simulations  is  not  optimized  for  state  estimation.  Still,  the  state  estimator  s  (TFiiter  &rid 
ajrue  values  are  less  than  their  corresponding  blended  values.  Both  horizontal  values  of 
Opiiter  converge  to  2.1  ft  in  the  beginning  of  the  scenario.  Then  they  gradually  increase  as 
a  increases  over  time.  There  is  a  slight  decrease  in  (JFUter  at  the  end,  reflecting  the 
decrease  in  interference  levels.  The  values  of  Ofrue  are  much  smaller  than  OFUter,  erring  on 
the  side  of  caution.  The  vertical  position  error  estimates  show  only  minor  improvement. 
This  can  be  attributed  to  the  radar  altimeter.  It  provides  better  measurements  as  the 
aircraft  approaches  the  ground  (viz.  Rnad),  making  both  the  MMAE’s  blended  output  and 
the  state  estimator’s  output  comparable  and  both  very  good. 

Tables  11-13  show  the  values  of  OFUter  and  arme  at  each  DH  for  the  three  position 
errors.  The  values  shown  are  greater  than  those  in  Case  I,  but  in  both  the  horizontal  and 
vertical  directions,  the  values  of  GFiuer  and  OTme  still  satisfy  the  Category  HI  requirements. 

Table  11.  Case  III:  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

OFOter  (ft) 

Ornie  (ft) 

I 

200 

3.567 

0.431 

n 

100 

3.868 

0.393 

m 

50 

3.841 

0.489 

Table  12.  Case  III:  North  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

OmUer  (ft) 

OTrue  (ft) 

I 

200 

3.672 

1.741 

n 

100 

3.976 

1.519 

III 

50 

3.966  n 

1.503 
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Table  13.  Case  III:  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

(^Filter  (ft) 

CTrue  (ft) 

I 

200 

1.902 

1.390 

II 

100 

1.166 

0.774 

m 

50 

0.714 

0.787 

The  level  of  interference  (measurement  noise)  was  significant  in  the  time  just  before  the 
DH  were  reach,  though  it  dropped  to  a  nominal  level  in  the  last  10  seconds  of  the 
scenario.  The  value  of  a  did  drop  dramatically  in  the  last  5  seconds,  so  this  may  have 
allowed  the  ILS  categories  to  be  met.  Case  V  shows  the  same  navigation  configuration 
under  a  continually  increasing  level  of  noise  to  observe  its  performance  when  the 
interference  level  is  high  at  touch-down. 

4.7.2.  Case  IV 

This  case  is  identical  to  Case  HI,  except  a  4.0  nmi/hr  drift  rate  INS  is  used  rather  than 
the  0.4  nmi/hr  INS.  The  results  of  d  and pk  from  the  MMAE  essentially  duplicate  those 
of  Case  m  shown  in  Figure  18,  as  expected,  thus  they  are  not  shown.  Figure  22  shows 
the  blended  output  of  the  position  errors  from  the  MMAE.  Like  the  results  of  Case  in 
shown  in  Figure  19,  these  plots  follow  the  results  of  White’s  research  [56]  (cf.  Nav  Case 
4  in  [56];  N.B.,  White  does  not  use  GPS  pseudolite  or  radar  altimeter  measurements  in 
this  case.).  Figure  23  shows  the  position  error  plots  from  the  output  of  the  ai  elemental 
filter  to  display  the  poor  estimates  generated  by  a  single  Kalman  filter  tuned  for  nominal 
performance. 
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Figure  22.  Case  IV:  Blended  (MMAE)  Position  Errors 


The  output  of  the  M^AE’s  state  estimator  is  shown  in  Figure  24.  As  in  Case  III,  there 
is  great  deal  of  improvement  of  these  position  errors  estimates  over  those  from  the 
MMAE,  especially  in  the  horizontal  directions.  The  values  of  Gfiuer  3nd  Ofrue  for  the 
horizontal  position  errors  are  significantly  greater  than  those  of  Case  III;  this  shows  the 
effects  of  using  an  INS  with  a  larger  drift  rate.  (Note  the  difference  in  scale  of  the  east 
and  north  position  error  axes  between  Figures  21  and  24.) 

There  is  a  noticeable  growth  in  these  values  after  3,850  seconds  into  the  scenario. 

The  is  expected  since  a  jumps  to  the  levels  of  800, 1600,  and  2000.  Also,  for  the  most 
part,  the  values  of  the  estimated  east  and  north  errors  ±  Grrue  still  fit  within  the  bounds  of 
(jputer-  The  vertical  position  error  estimates  do  not  show  significant  difference  between 
the  state  estimator  output  of  Case  HI,  though,  as  mentioned  previously,  the  radar  altimeter 
measurements  have  a  significant  impact  on  the  altitude  measurements. 

Tables  14  -  16  show  the  values  of  Gputer  and  aime  at  each  DH  for  the  three  position 
errors.  In  comparing  the  values  of  Gputer  against  the  FAA  criteria,  this  case  still  satisfies 
Category  IE  for  the  vertical  positioning  error;  however,  it  can  only  meet  Category  I  for 
horizontal  positioning.  The  values  of  Gprue  show  better  results  with  vertical  positioning 
meeting  Category  HI  and  horizontal  positioning  meeting  Category  11.  Unlike  Case  HI, 
the  horizontal  position  errors  did  not  improve  quickly  enough  to  met  Category  III  in  all 
areas.  It  can  be  seen  that  the  system  modeled  with  the  4.0  nmi/hr  INS  results  in  degraded 
performance  in  horizontal  position.  Thus,  this  suggests  that  an  INS  with  a  lower  drift  rate 
is  important  to  satisfy  the  landing  requirements  in  this  scenario. 
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Figure  24.  Case  IV:  State  Estimator  Position  Errors 


90 


Table  14.  Case  IV:  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

Gtme  (ft) 

I 

200 

10.67 

5.614 

II 

100 

11.72 

5.579 

III 

50 

10.43 

3.394 

Table  15.  Case  IV:  North  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

OTrue  (ft) 

I 

200 

11.42 

8.222 

n 

100 

12.58 

5.980 

m 

50 

11.25 

4.033 

Table  16.  Case  IV:  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

CFrrue  (ft) 

I 

200 

2.016 

1.954 

n 

100 

1.273 

0.920 

m 

50 

0.798 

0.818 

4.8.  Increasing  Noise  Levels  Scenarios 

Cases  V  and  VI  examine  the  system  performance  against  a  scenario  where  the  level 
of  GPS  pseudorange  measurement  noise  increases  as  the  aircraft  approaches  the  touch¬ 
down  point.  The  level  of  a  increases  in  inverse  proportion  to  the  square  of  the  distance  to 
the  RF  interference  source  (viz.  Friis  transmission  formula  [8]).  The  source  is  assumed  to 
be  in-line  with  the  flight  path  of  the  aircraft.  This  is  meant  to  depict  a  point  source  of  RF 
interference  increasing  the  level  of  measurement  noise  for  Rpuj.  Although  this  is  not  the 
most  accurate  representation,  it  is  meant  to  show  the  effects  on  the  position  errors  as  the 
landing  aircraft  moves  toward  an  interference  source. 
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To  implement  the  increasing  level  of  3  over  time  (hence  over  distance,  as  the  aircraft 
maintains  a  steady  glide  velocity),  different  constant  values  of  3  had  to  be  shown  in 
separate  intervals.  This  is  due  the  limitation  of  the  software  configuration  used,  which 
could  only  accept  constant  parameter  values  at  a  maximum  of  19  intervals.  Table  17 
shows  the  how  3  varies  over  this  scenario. 


Table  17.  Increasing  Noise  Levels  Scenario 


Interval 

Actual 

Start  Time 

Parameter 

(sec) 

(a) 

3,810 

1.0 

3,820 

8.0 

3,825 

30.0 

3,830 

70.0 

3,835 

125.0 

3,840 

195.0 

3,845 

280.0 

3,850 

380.0 

3,855 

500.0 

3,860  n 

630.0 

3,865 

780.0 

3,870 

945.0 

3,875 

1,125.0 

3,880 

1,320.0 

3,885 

1,530.0 

3,890 

1,760.0 

3,895 

2,000.0 

4.8.1.  Case  V 

This  case  involves  the  navigation  system  using  the  0.4  nmi/hr  drift  rate  INS.  Figure 
25  shows  how  d  and  pk  vary  over  the  scenario.  As  in  Figure  18,  the  top  plot  of  Figure  25 
includes  the  actual  value  of  the  parameter  a  with  heavy  solid  lines.  The  value  of  d 
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follows  a  throughout  the  scenario.  Up  to  3,820  seconds,  the  a\  elemental  filter  is  chosen 
to  be  the  best  model.  When  a  =  8,  between  3,820  seconds  and  3,825,  the  MMAE 
calculates  a  blended  estimate  of  d  =  93.  After  that  point,  when  a  becomes  30  at  3,825 
seconds,  the  MMAE  gives  nearly  all  of  the  probability  weighting  to  the  a2  elemental 
filter.  This  remains  the  case  until  3,850  seconds,  when  a  becomes  380.  From  this  point 
on,  blending  occurs  within  the  MMAE  to  determine  the  value  d.  At  3,875  seconds,  the 
estimated  parameter,  d,  starts  to  lag  the  actual  parameter,  a,  which  takes  on  the  value  of 
1,125  at  that  time.  Throughout  the  remainder  of  the  scenario,  the  value  of  d  slowly 
increases  towards  its  upper  limit,  just  under  as.  Again,  this  highlights  the  impact  of  not 
having  an  elemental  filter  at  a  =  2,000  (see  recommendation  in  Section  5.3). 

Figure  26  shows  the  position  errors  according  to  the  M^AE’s  state  estimator.  The 
results  appear  similar  to  those  of  Case  m  shown  in  Figure  21  without  the  decrease  in 
Gf  liter  at  the  very  end.  The  absence  of  the  decrease  is  due  to  a  remaining  at  2,000,  rather 
taking  on  lower  levels  as  in  Case  HI.  Also,  the  values  of  OTme  are  much  less  than  Gputer  as 
in  Case  HI,  though  they  are  greater  than  Ufrue  of  Case  m  overall. 

Tables  18-20  show  the  values  of  Gfiuer  and  OTrue  at  each  DH.  The  horizontal  and 
vertical  positioning  error  according  to  both  Gpiuer  and  Grrue  continues  to  meet  Category 
III.  This  suggests  there  is  no  significant  impact  of  a  remaining  at  2,000  in  this  scenario. 
Thus,  the  navigation  system  configuration  with  the  lower  drift  rate  INS  can  operate 
sufficiently  well  in  the  presence  of  a  high  degree  of  pseudorange  measurement  noise. 
However,  it  should  noted  that  GPS  aiding  was  present  at  the  beginning  of  the  scenario 
and  reduced  over  time.  Previous  research  [7, 13,  38,  56]  has  shown  that,  without  GPS 
aiding,  this  level  of  precision  could  not  be  met. 
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Table  18.  Case  V:  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

Opilter  (ft) 

(TTrue  (ft) 

I 

200 

3.566 

0.863 

II 

100 

3.914^ 

0.983 

III 

50  1 

4.153 

1.149 

Table  19.  Case  V;  North  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

(^Filter  (ft) 

OTrue  (ft) 

I 

200 

3.649 

1.481 

II 

100 

3.997 

1.309 

m 

50  1 

4.236  ^ 

1.712 

Table  20.  Case  V:  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

OTrue  (ft) 

I 

200 

1.902 

1.395 

n 

100 

1.166 

0.768 

m 

50  1 

0.715 

0.794 

4.8.2.  Case  VI 

This  case  is  identical  to  Case  V,  except  for  the  4.0  nmi/hr  drift  rate  INS.  As  in  Case 
rv,  the  results  of  d  and  pk  are  not  shown  since  they  essentially  duplicate  those  of  the 
previous  case  shown  in  Figure  25.  The  output  of  the  state  estimator  is  shown  in  Figure 
27.  The  standard  deviations  of  the  horizontal  position  errors  are  significantly  greater  than 
those  of  Case  V.  These  are  in  the  same  scale  as  in  Case  IV,  so  Figure  27  uses  the  same 
axis  scaling  as  Figure  24.  Unlike  Case  IV,  the  values  of  CpiUer  and  Oprue  for  the  horizontal 
errors  continue  to  increase  due  to  the  growth  in  the  level  of  a,  in  the  same  manner  as  seen 
in  Case  V. 
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Tables  21-23  show  the  values  of  Gputer  and  aprue  the  DH  points.  As  in  Case  IV,  the 
value  of  Ofiiter  meets  Category  III  in  the  vertical  direction  and  only  Category  I  in  the 
horizontal  directions.  The  values  of  Oprue  meet  Category  III  specifications  in  the  vertical 
direction;  however,  they  can  only  meet  Category  II  specifications  in  the  horizontal 
directions.  The  filter-computed  standard  deviation,  Gputer,  is  more  conservative  than  the 
standard  deviation  over  the  Monte  Carlo  runs,  Gprue,  as  usual.  As  in  Case  IV,  the  quality 
of  the  INS  has  a  prominent  effect  on  the  position  errors  in  the  presence  of  measurement 
noise.  This  also  emphasizes  the  importance  of  having  an  INS  of  good  quality  on  board 
when  GPS  measurements  are  disrupted. 

Table  21.  Case  VI:  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

GpiUer  (ft) 

GTrue  (ft) 

I 

200 

14.46 

7.724 

n  ^ 

100 

15.86 

6.208 

m 

50  1 

16.68  ^ 

8.983 

Table  22.  Case  VI:  North  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

GpiUer  (ft) 

Gprue  (ft) 

I 

200 

15.23 

10.29 

n 

100 

16.81 

7.921 

m 

50  n 

17.74 

12.64 

Table  23.  Case  VI:  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

(^Filter  (ft) 

G'True  (ft) 

I 

200 

2.017 

1.942 

n  n 

100 

1.213 

0.923 

m 

50 

0.800  n 

0.829 
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4.9.  Decreasing  Noise  Levels  Scenarios 

In  the  last  two  cases,  the  system  performance  is  examined  against  a  scenario  where 
the  level  of  GPS  pseudorange  measurement  noise  starts  at  the  level  where  GPS  is 
ineffective.  As  the  trajectory  progresses,  the  interference  decreases  to  a  nominal  level  as 
the  aircraft  approaches  the  touch-down  point.  The  level  of  a  decreases  in  inverse 
proportion  to  the  square  of  the  distance  to  the  RF  interference  source,  as  in  Cases  V  and 
VI.  Also,  the  source  is  assumed  to  be  in-line  with  the  flight  path  of  the  aircraft.  Table  24 
shows  the  how  a  varies  over  this  scenario. 


Table  24.  Decreasing  Noise  Levels  Scenario 


Interval 

Actual 

Start  Time 

Parameter 

(sec) 

(a) 

3,810 

2,000.0 

3,820 

1,760.0 

3,825 

1,530.0 

3,830 

1,320.0 

3,835 

1,125.0 

3,840 

945.0 

3,845 

780.0 

3,850 

630.0 

3,855 

500.0 

3,860 

380.0 

3,865 

280.0 

3,870 

195.0 

3,875 

125.0 

3,880 

70.0 

3,885 

30.0 

3,890 

8.0 

3,895 

1.0 
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4.9.1.  Case  VII 

As  with  the  other  odd-numbered  cases,  this  case  involves  the  navigation  system  using 
the  0.4  nmi/hr  drift  rate  INS.  Figure  28  shows  how  d  and  pk  vary  over  the  scenano.  As 
in  Case  V,  the  value  of  a  follows  a  throughout  the  scenario  with  a  lag.  The  effect  of  this 
is  that  d  tends  to  be  greater  than  a,  except  during  the  beginning  of  the  scenario  (d  tended 
to  be  less  than  a  in  previous  scenario  where  a  was  increasing  due  to  as  <  2,000).  The 
value  of  d  reaches  the  upper  bound  just  under  as  by  the  time  a  drops  to  1,530.  Up  to  the 
time  of  3,880  seconds  in  the  scenario,  blending  is  achieved  within  the  MMAE  as  the 
probability  plots  of  the  az,  as,  a^,  and  as  elemental  filters  show.  From  3,880  to  3,895 
seconds  of  the  scenario,  the  MMAE  favors  the  az  elemental  filter  almost  exclusively,  as 
in  Case  V,  when  the  values  of  a  was  between  ai  and  az.  After  that  point,  when  a  reaches 
the  nominal  level  of  1,  the  MMAE  readjusts  d  within  5  seconds. 

Figure  29  shows  the  positions  errors  according  to  the  M^AE’s  state  estimator.  The 
values  of  apuur  and  CTirue  decrease,  as  the  state  estimator  adapts  to  account  for  the 
decrease  in  a  as  the  scenario  progresses.  Note  that  the  scale  of  the  horizontal  error  plots 
is  greater  than  in  previous  figures,  to  show  the  entire  range  of  how  the  errors  and  their 
standard  deviations  vary.  Overall,  the  values  of  Crrue  stay  within  the  aniter  bounds, 
showing  that  the  state  estimator  filter  performs  adequately. 

Tables  25  -  27  show  the  values  of  OFUter  and  anue  at  each  DH.  For  both  the  Cpiuer  and 

OTrue  in  vertical  and  horizontal  directions.  Category  HI  is  satisfied.  Despite  a  high  level  of 
interference  at  the  beginning  of  the  scenario,  the  state  estimator,  provided  with  updated 
values  of  d,  can  perform  well  enough  to  meet  Category  HI.  This  indicates  that  the  effects 
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Parameter  Value 


Vertical  Error  (ft)  North  Error  (ft)  East  Error  (ft) 


State  Estimator  Output 


Time  (sec) 


Figure  29.  Case  VII:  State  Estimator  Position  Errors 
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Table  25.  Case  VII:  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

^True  (ft) 

I 

200 

5.507 

0.936 

n 

100 

1.858 

0.701 

m 

50 

1.509 

0.582 

Table  26.  Case  VII:  North  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

OpiUer  (ft) 

Orrue  (ft) 

I 

200 

5.982 

3.183 

n 

100 

2.068 

0.888 

m 

50 

1.665 

0.838 

Table  27.  Case  VII:  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

^Filter  (ft) 

(TTrue  (ft) 

I 

200 

1.886 

1.435 

n 

100 

1.146 

0.719 

m 

50 

0.710 

0.798 

of  increased  pseudorange  measurement  noise  do  not  linger  within  the  filters’  calculations 
long  enough  to  affect  the  outcomes,  if  the  noise  variance  returns  to  a  nominal  level  soon 
before  touch-down. 

4.9.2.  Case  VIII 

This  case  is  identical  to  Case  VII,  except  for  the  4.0  nmi/hr  drift  rate  INS.  As  with 
the  other  even-numbered  cases,  the  plots  of  a  and  pk  are  essentially  identical  to  its 
previous  case  (see  Figure  28).  Figure  30  shows  the  positions  errors  determined  by  the 
M^AE’s  state  estimator.  The  only  notable  feature  is  that  OFUter  for  the  horizontal  position 
errors  remains  at  a  higher  level  than  in  Case  VII  from  3,870  to  3,890  seconds  in  the 


103 


Vertical  Error  (ft)  North  Error  (ft)  East  Error  (ft) 


State  Estimator  Output 


Time  (sec) 


Figure  30.  Case  VIII:  State  Estimator  Position  Errors 
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scenario.  Otherwise,  the  navigation  system  demonstrates  a  performance  similar  to  that  in 
Case  VII. 

Tables  28  -  30  show  the  values  of  (ypnter  and  axme  at  the  DH  points.  Although  the 
values  of  OFUter  and  (JTme  at  the  DH  points  are  greater  than  in  Case  VII,  Category  III  is 
satisfied,  in  both  the  vertical  and  horizontal  directions.  As  with  Case  VII,  the  effects  of 
increased  measurement  noise  do  not  persist  if  a  nominal  level  is  reached  in  reasonably 
short  time  prior  to  touch-down. 

Table  28.  Case  VIII;  East  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

(^Filter  (ft) 

CTlrue  (ft) 

I 

200 

8.965 

2.495 

n 

100 

2.691 

m 

50 

2.210 

■if.-ltiiM 

Table  29.  Case  VIII;  North  Position  Error  Standard  Deviations 


ILS  Category 

OFUter  (ft) 

^True  (ft) 

I 

WKESSam 

2.713 

II 

100 

2.979 

1.767 

m 

50 

2.423 

1.524 

Table  30.  Case  VIII;  Vertical  Position  Error  Standard  Deviations 


ILS  Category 

DH  (ft) 

(^Filter  (ft) 

(^True  (ft) 

I 

200 

n 

100 

MEm 

III 

50 

0.793 

0.825 
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4.10.  Summary 

The  test  cases  presented  in  this  chapter  covered  four  different  pseudorange 
measurement  noise  scenarios  using  the  flight  profile  of  a  simulated  landing  for  a  KC-135 
tanker  aircraft.  Table  31  summarizes  the  results  of  all  the  test  cases.  The  minimum  ILS 
category  each  case  can  satisfy  is  listed,  i.e.,  the  worst  performance  between  the  horizontal 
and  vertical  position  errors,  and  between  the  filter-computed  standard  deviation,  opiUeri 
and  the  standard  deviation  over  the  Monte  Carlo  runs,  Crme,  is  used  to  make  this 
assessment.  For  most  instances.  Category  HI  could  be  met.  In  the  cases  which  could  not, 
the  higher  drift  rate  INS  were  used  and  there  was  a  high  level  of  interference  near  the  DH 
points. 


Table  31.  ILS  Precision  Approach  Categories  Met  for  each  Test  Case 


ILS 

Case 

Case 

Case 

Case 

Case 

Case 

Case 

Case 

Category 

I 

II 

III 

IV 

V 

VI 

VII 

VIII 

I 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

n 

Yes 

Yes 

Yes 

No 

Yes 

No 

Yes 

Yes 

m 

Yes 

Yes 

Yes 

No 

Yes 

No 

Yes 

Yes 
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5.  Conclusions  and  Recommendations 


5.1.  Overview 

The  objective  of  this  research  effort  has  been  to  apply  the  M^AE  architecture  to  an 
aircraft  precision  approach  landing  where  the  GPS  signals  were  subjected  to  interference. 
The  work  has  focused  on  calculating  accurate  data  on  the  aircraft’s  position  throughout 
its  glideslope,  and  comparing  the  precision  of  the  data  against  FAA ILS  criteria  to  assess 
the  performance  of  the  modeled  aircraft  navigation  system.  The  navigation  system 
models  have  included  an  INS  with  a  built-in  barometric  altimeter,  a  GPS  receiver 
collecting  measurements  from  4  GPS  SVs  and  a  ground-based  pseudolite,  and  a  radar 
altimeter.  This  is  a  direct  extension  to  the  research  of  previous  AFTT  theses  by  Gray  [13], 
Britton  [7],  White  [56],  and  Miller  [33]. 

The  modeled  system  underwent  a  series  of  simulations  where  the  measurement  noise 
of  the  GPS  signals  varied  to  depict  RF  interference.  An  M^AE  architecture  was  applied 
to  handle  the  nature  of  the  changing  parameter  of  the  noise,  yet  still  produce  accurate 
state  estimation.  This  has  been  done  to  attempt  to  maintain  the  level  of  precision  in  the 
estimated  aircraft’s  position. 

5.2.  Conclusions 

This  research  showed  that  an  M^AE  provides  improvement  in  estimating  both  GPS 
interference  levels  and  aircraft  position  errors  over  that  obtained  in  previous  research. 
Several  test  cases  were  researched  using  two  INSs  of  different  quality  and  various 
interference  profiles.  The  overall  performance  of  each  simulation  was  judged  against  the 
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FAA ILS  criteria  for  precision  approach.  Most  of  the  test  cases  examined  in  this  research 
effort  satisfied  the  requirements  for  a  Category  III  approach.  In  the  cases  where  the 
Category  III  requirements  were  not  met,  there  was  a  high  level  of  GPS  interference  near 
the  landing  point  and  the  navigation  system  was  using  the  poorer  quality  INS  (4.0  nmi/hr 
drift  rate).  However,  in  those  cases.  Category  I  requirements  were  met.  When  Category 
m  could  not  be  met,  it  was  due  to  the  level  of  accuracy  in  the  horizontal  (east  and  north) 
directions.  In  the  horizontal  directions,  the  INS  has  aiding  from  GPS  measurements  only, 
unlike  in  the  vertical  direction,  where  there  are  measurements  from  a  barometric  altimeter 
and  a  radar  altimeter  as  well  GPS  (recall  that  INSs  inherently  have  a  greater  difficulty 
estimating  vertical  positions  [6, 27],  an  extremely  important  concern  to  all  aircraft,  thus 
these  instruments  are  added).  The  results  of  the  simulations  show  the  significance  of  the 
INS  drift  rate  and  GPS  measurements.  Without  high-accuracy  pseudorange 
measurements  provided  by  the  GPS  receiver,  the  effect  of  the  INS  drift  rate  becomes 
prominent  enough  to  affect  the  positioning  accuracy  at  the  decision  height  (DH)  points. 

Based  on  the  results,  the  following  conclusions  can  be  made: 

1.  The  M^AE’s  MMAE  provided  good  estimates  of  the  interference  levels,  except  at 
the  high  interference  levels  (beyond  the  boundary  of  the  MMAE’s  elemental 
filters)  where  parameter  estimation  suffered.  This  would  be  readily  compensated 
by  allowing  an  elemental  filter  to  assume  the  largest  interference  level  of  concern. 
However,  navigation  performance  was  not  adversely  affected. 

2.  The  M^AE  architecture  and  techniques  were  validated  against  a  realistic  aircraft 
navigation  truth  model.  The  results  in  this  research  provide  a  practical 
demonstration  of  M^AE  using  truth  and  filter  models  of  different  orders. 
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3.  The  level  of  GPS  pseudorange  measurement  noise  had  a  greater  impact  on  the 
position  errors  in  the  horizontal  directions  (east  and  north)  than  in  the  vertical 
direction.  This  is  attributed  to  the  assistance  the  navigation  system  derives  from 
other  altitude  measuring  devices  that  are  unaffected  by  interference  in  the  GPS 
frequency  band,  particularly  the  very  accurate  radar  altimeter. 

5.3.  Recommendations 

The  scope  of  this  research  encompassed  examining  the  performance  of  a  GPS-aided 
INS  in  the  presence  of  GPS  interference.  As  shown,  precision  approach  under  the 
presence  of  GPS  interference  was  achieved  using  an  M^AE  architecture.  However,  there 
are  several  areas  that  require  further  research: 

1.  Develop  a  more  accurate  model  of  interference  than  the  one  implemented  here. 
The  effects  of  a  change  in  the  actual  power  of  RF  noise  should  be  represented  by 
a  more  detailed  algorithm. 

2.  Implement  a  moving-bank  MMAE  architecture,  developed  by  Vasquez  [54,  55], 
in  place  of  the  fixed-bank  MMAE  parameter  estimator  to  investigate  the 
differences  in  performance. 

3.  Generate  several  different  flight  profiles  to  examine  the  effects  of  maneuvering 
during  the  landing  approach  phase  of  flight. 

4.  Perform  the  constrained  parameter  discretization  fixing  filters  at  ai  =  1  and  = 
2,000  (where  there  are  K  filters  in  the  MMAE). 

5.  Incorporate  logic  where,  if  =  1  -  (K-l)pmin,  then  pk=landpj  =  0y  j*k3]e 
{ 1,  2, ...,  AT}  is  sent  the  output  of  the  MMAE. 
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Appendix  A.  Acronym  List 


AFTT 

Air  Force  Institute  of  Technology 

AGL 

Above  Ground  Level 

DGPS 

Differential  GPS 

DH 

Decision  Height 

ECEF 

Earth-Centered,  Earth-Fixed  (reference  frame) 

EKE 

Extended  Kalman  Filter 

FAA 

Federal  Aviation  Administration 

GPS 

Global  Positioning  System 

ILS 

Instrument  Landing  System 

INS 

Inertial  Navigation  System 

IRDF 

Inter-Residual  Distance  Feedback 

LAAS 

Local  Area  Augmentation  System 

MSOFE 

Multi-mode  Simulation  of  Optimal  Filter  Equations 

MMSOFE 

Multiple  Model  Simulation  of  Optimal  Filter  Equations 

MMAE 

Multiple  Model  Adaptive  Estimator 

M^AE 

Modified  Multiple  Model  Adaptive  Estimator 

PLS 

Precision  Landing  System 

PPS 

Precise  Positioning  Service 

PR 

Pseudorange 

RF 

Radio  Frequency 

SA 

Selective  Availability 

SP 

Standard  Positioning  Service 

SV 

Space  Vehicle 

TACAN 

Tactical  Air  Navigation 

UHF 

Ultra  High  Frequency 

VHF 

Very  High  Frequency 

VOR 

VHF  Omnidirectional  Range 

WGN 

White  Gaussian  Noise 

WGS-84 

World  Geodetic  Survey  1984 
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Appendix  B.  Error  States  of  the  Truth  and  Filter  Models 


This  appendix  lists  the  error  states  implemented  by  the  truth  and  filter  models 
considered  in  this  research.  Tables  32  -  34  show  a  complete  list  of  the  93  error  states  that 
comprise  the  Litton  LN-93  BSTS  truth  model  [26].  Tables  35  and  36  list  the  39  error  states 
of  the  reduced-order  truth  model  used  in  the  simulations  [7, 13, 38,  56,  57].  These 
consist  of  a  proper  subset  of  the  original  LN-93  states,  i.e.,  every  state  in  the  reduced- 
order  truth  model  is  a  state  in  the  original  truth  model.  Table  37  lists  the  30  error  states 
characterizing  a  GPS  truth  model  using  4  SVs  and  1  receiver  [7, 13,  38,  56,  57].  Table 
38  lists  the  22  error  states  associated  with  a  DGPS  truth  model;  these  form  a  proper 
subset  of  the  GPS  states  [7, 13,  38, 56, 57].  Table  39  lists  the  13-state  filter  model, 
consisting  of  11  INS  error  states  and  the  2  error  states  related  to  the  GPS  receiver  [7, 13, 
38, 56,  57]. 


Ill 


Table  32.  Original  93-State  Truth  Model  for  the  INS:  States  1-31 


State 

1^9 

State 

Number 

DeHnition 

1 

5e, 

X-component  of  vector  angle  from  true  to  computer  frame 

2 

SOy 

Y-component  of  vector  angle  from  true  to  computer  frame 

3 

5d, 

Z-component  of  vector  angle  from  true  to  computer  frame 

4 

X-component  of  vector  angle  from  true  to  platform  frame 

5 

(th 

Y-component  of  vector  angle  from  true  to  platform  frame 

6 

<1>Z 

Z-component  of  vector  angle  from  true  to  platform  frame 

7 

X-component  of  error  in  computed  velocity 

8 

dVy 

Y-component  of  error  in  computed  velocity 

9 

sv. 

Z-component  of  error  in  computed  velocity 

10 

dh 

Error  in  vehicle  altitude  above  reference  ellipsoid 

11 

mm 

Error  in  lagged  inertial  altitude 

12 

SSs 

Error  in  vertical  channel  aiding  state 

13 

SS4 

Error  in  vertical  channel  aiding  state 

14 

K 

X-component  of  gyro  correlated  drift  rate 

15 

K 

Y-component  of  gyro  correlated  drift  rate 

K 

Z-component  of  gyro  correlated  drift  rate 

17 

V. 

X-component  of  accelerometer  and  velocity  quantizer  correlated  noise 

18 

Y-component  of  accelerometer  and  velocity  quantizer  correlated  noise 

19 

Z-component  of  accelerometer  and  velocity  quantizer  correlated  noise 

■B91I 

8gx 

X-component  of  gravity  vector  errors 

21 

Sgy 

Y-component  of  gravity  vector  errors 

22 

Z-component  of  gravity  vector  errors 

23 

ShB 

Total  barometric  altimeter  correlated  error 

24 

K 

X-component  of  gyro  trend 

25 

K 

Y-component  of  gyro  trend 

26 

K 

Z-component  of  gyro  trend 

27 

V. 

X-component  of  accelerometer  trend 

28 

Y-component  of  accelerometer  trend 

29 

Z-component  of  accelerometer  trend 

30 

bx 

X-component  of  gyro  drift  rate  repeatability 

31 

_ h _ 

Y-component  of  gyro  drift  rate  repeatability 

Table  33.  Original  93-State  Truth  Model  for  the  INS:  States  32  -  63 


State 

Number 

State 

Definition 

32 

Z-component  of  gyro  drift  rate  repeatability 

33 

X-component  of  gyro  scale  factor  error 

34 

Y-component  of  gyro  scale  factor  error 

35 

Sg? 

Z-component  of  gyro  scale  factor  error 

36 

mm 

X  gyro  misalignment  about  Y-axis 

37 

Y  gyro  misalignment  about  X-axis 

38 

X3 

Z  gyro  misalignment  about  X-axis 

39 

V] 

X  gyro  misalignment  about  Z-axis 

V2 

Y  gyro  misalignment  about  Z-axis 

41 

V3 

Z  gyro  misalignment  about  Y-axis 

42 

Dxxx 

X  gyro  scale  factor  nonlinearity 

43 

Y  gyro  scale  factor  nonlinearity 

44 

Z  gyro  scale  factor  nonlinearity 

45 

X  gyro  scale  factor  asymmetry  error 

46 

— 

Y  gyro  scale  factor  asymmetry  error 

47 

■■ 

Z  gyro  scale  factor  asymmetry  error 

48 

X-component  of  accelerometer  bias  repeatability 

49 

Y-component  of  accelerometer  bias  repeatability 

Z-component  of  accelerometer  bias  repeatability 

51 

Sax 

X-component  of  accelerometer  and  velocity  quantizer  scale  factor 

52 

Y-component  of  accelerometer  and  velocity  quantizer  scale  factor 

53 

Saz 

Z-component  of  accelerometer  and  velocity  quantizer  scale  factor 

54 

Sqax 

X-component  of  accelerometer  and  velocity  quantizer  scale  factor  asymmetry 

55 

Y-component  of  accelerometer  and  velocity  quantizer  scale  factor  asymmetry 

56 

SoAz 

Z-component  of  accelerometer  and  velocity  quantizer  scale  factor  asymmetry 

57 

fxx 

Coefficient  of  error  proportional  to  square  of  measured  acceleration 

58 

Coefficient  of  error  proportional  to  square  of  measured  acceleration 

59 

Coefficient  of  error  proportional  to  square  of  measured  acceleration 

Coefficient  of  error  proportional  to  square  of  measured  acceleration 

61 

Coefficient  of  error  proportional  to  square  of  measured  acceleration 

62 

iB& 

Coefficient  of  error  proportional  to  square  of  measured  acceleration 

63 

imsm 

Coefficient  of  error  proportional  to  square  of  measured  acceleration 

Table  34.  Original  93-State  Truth  Model  for  the  INS:  States  64-93 


State 

Number 

H^IQI 

State 

Deflnition 

64 

fzz 

Coefficient  of  error  proportional  to  square  of  measured 
acceleration 

65 

fzy 

Coefficient  of  error  proportional  to  square  of  measured 
acceleration 

66 

til 

X  accelerometer  misalignment  about  Z-axis 

67 

ti2 

Y  accelerometer  misalignment  about  Z-axis 

68 

tl3 

Z  accelerometer  misalignment  about  Y-axis 

69 

03 

Z  accelerometer  misalignment  about  Y-axis 

70 

X-component  of  accelerometer  bias  thermal  transient 

71 

Y-component  of  accelerometer  bias  thermal  transient 

72 

Z-component  of  accelerometer  bias  thermal  transient 

73 

K 

X-component  of  initial  gyro  drift  rate  bias  thermal  transient 

74 

K 

Y-component  of  initial  gyro  drift  rate  bias  thermal  transient 

75 

K 

Z-component  of  initial  gyro  drift  rate  bias  thermal  transient 

76 

X  gyro  compliance  term 

77 

lESH 

X  gyro  compliance  term 

78 

— 

X  gyro  compliance  term 

79 

X  gyro  compliance  term 

80 

P xzz 

X  gyro  compliance  term 

81 

Fxz 

X  gyro  compliance  term 

82 

WKsKm 

Y  gyro  compliance  term 

83 

Y  gyro  compliance  term 

84 

Y  gyro  compliance  term 

85 

wtmm 

Y  gyro  compliance  term 

86 

Y  gyro  compliance  term 

87 

Y  gyro  compliance  term 

88 

mmm 

Z  gyro  compliance  term 

89 

Z  gyro  compliance  term 

90 

BE1 

Z  gyro  compliance  term 

91 

Z  gyro  compliance  term 

92 

Z  gyro  compliance  term 

93 

HSflli 

Z  gyro  compliance  term 

Table  35.  Reduced-Order  39-State  Truth  Model  for  the  INS;  States  1-20 


State 

Number 

l^H 

State 

DeHnition 

1 

X-component  of  vector  angle  from  true  to  computer  frame 

2 

59, 

Y-component  of  vector  angle  from  true  to  computer  frame 

3 

59, 

Z-component  of  vector  angle  from  true  to  computer  frame 

4 

4 

X-component  of  vector  angle  from  true  to  platform  frame 

5 

Y-component  of  vector  angle  from  true  to  platform  frame 

6 

Z-component  of  vector  angle  from  true  to  platform  frame 

7 

— 

X-component  of  error  in  computed  velocity 

8 

Y-component  of  error  in  computed  velocity 

9 

Z-component  of  error  in  computed  velocity 

10 

Error  in  vehicle  altitude  above  reference  ellipsoid 

11 

Shs 

Total  barometric  altimeter  correlated  error 

12 

Error  in  lagged  inertial  altitude 

13 

Error  in  vertical  channel  aiding  state 

14 

Error  in  vertical  channel  aiding  state 

15 

V. 

X-component  of  accelerometer  and  velocity  quantizer  correlated  noise 

16 

Y-component  of  accelerometer  and  velocity  quantizer  correlated  noise 

17 

Z-component  of  accelerometer  and  velocity  quantizer  correlated  noise 

18 

Sgx 

X-component  of  gravity  vector  errors 

19 

Y-component  of  gravity  vector  errors 

5^. 

Z-component  of  gravity  vector  errors 
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Table  36.  Reduced-Order  39-State  Truth  Model  for  the  INS:  States  21  -  39 


State 

1^9 

State 

Number 

DeHnition 

21 

bx 

X-component  of  gyro  drift  rate  repeatability 

22 

by 

Y-component  of  gyro  drift  rate  repeatability 

23 

bz 

Z-component  of  gyro  drift  rate  repeatability 

24 

X-component  of  gyro  scale  factor  error 

25 

Y-component  of  gyro  scale  factor  error 

26 

Z-component  of  gyro  scale  factor  error 

27 

X-component  of  accelerometer  bias  repeatability 

28 

Y-component  of  accelerometer  bias  repeatability 

29 

Z-component  of  accelerometer  bias  repeatability 

30 

Sax 

X-component  of  accelerometer  and  velocity  quantizer  scale  factor 

31 

KM 

Y-component  of  accelerometer  and  velocity  quantizer  scale  factor 

32 

Saz 

Z-component  of  accelerometer  and  velocity  quantizer  scale  factor 

33 

SoAx 

X-component  of  accelerometer  and  velocity  quantizer  scale  factor  asymmetry 

34 

liBiSilli 

Y-component  of  accelerometer  and  velocity  quantizer  scale  factor  asymmetry 

35 

KM 

Z-component  of  accelerometer  and  velocity  quantizer  scale  factor  asymmetry 

36 

X  accelerometer  misalignment  about  Z-axis 

37 

IH2 

Y  accelerometer  misalignment  about  Z-axis 

38 

H3 

Z  accelerometer  misalignment  about  Y-axis 

39 

(73 

Z  accelerometer  misalignment  about  X-axis 
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Table  37.  30-State  Truth  Model  for  GPS 


State 

State 

Number 

Definition 

1 

User  clock  bias 

2 

^^Uclk 

User  clock  drift 

3 

SV  1  multi-path  error 

4 

SV  1  tropospheric  error 

5 

SV  1  ionospheric  error 

6 

SV  1  clock  error 

7 

SV  1  x-component  of  position  error 

8 

Sysvi 

SV  1  y-component  of  position  error 

9 

^^SVl 

SV  1  z-component  of  position  error 

^^MP2 

SV  2  multi-path  error 

11 

SV  2  tropospheric  error 

12 

SR„2 

SV  2  ionospheric  error 

13 

SV  2  clock  error 

14 

Sx^vi 

SV  2  x-component  of  position  error 

15 

^SV2 

SV  2  y-component  of  position  error 

16 

^^SV2 

SV  2  z-component  of  position  error 

17 

^^AfP3 

SV  3  multi-path  error 

18 

SV  3  tropospheric  error 

19 

SV  3  ionospheric  error 

20 

SV  3  clock  error 

21 

^*^5V3 

SV  3  x-component  of  position  error 

22 

^5V3 

SV  3  y-component  of  position  error 

23 

^^SV3 

SV  3  z-component  of  position  error 

24 

SV  4  multi-path  error 

25 

^^tropA 

SV  4  tropospheric  error 

26 

^Pioni 

SV  4  ionospheric  error 

27 

III^SRII^I 

SV  4  clock  error 

28 

5^5^  4 

SV  4  x-component  of  position  error 

29 

^SV4 

SV  4  y-component  of  position  error 

30 

&5V4 

SV  4  z-component  of  position  error 
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Table  38.  22-State  Truth  Model  for  Differential  GPS 


State 

ISB 

State 

Number 

Definition 

1 

User  clock  bias 

2 

User  clock  drift 

3 

SV  1  tropospheric  error 

4 

SV  1  ionospheric  error 

5 

SV  1  x-component  of  position  error 

6 

S  V  1  y-component  of  position  error 

7 

SV  1  z-component  of  position  error 

8 

SR.r„2 

SV  2  tropospheric  error 

9 

SR^2 

SV  2  ionospheric  error 

10 

5x^y2 

SV  2  x-component  of  position  error 

11 

^SVl 

SV  2  y-component  of  position  error 

12 

&5y2 

SV  2  z-component  of  position  error 

13 

SV  3  tropospheric  error 

14 

SV  3  ionospheric  error 

15 

^'^5V3 

SV  3  x-component  of  position  error 

16 

SV  3  y-component  of  position  error 

17 

SV  3  z-component  of  position  error 

18 

^Kop* 

SV  4  tropospheric  error 

19 

SR^, 

SV  4  ionospheric  error 

20 

5Xsy4 

SV  4  x-component  of  position  error 

21 

^SV4 

SV  4  y-component  of  position  error 

22 

&SV4 

SV  4  z-component  of  position  error 
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Table  39.  13-State  Filter  Model  (for  all  navigation  components) 


State 

IEIS3I 

State 

Number 

Deflnition 

1 

56, 

X-component  of  vector  angle  from  true  to  computer  frame 

2 

Y-component  of  vector  angle  from  true  to  computer  frame 

3 

56, 

Z-component  of  vector  angle  from  true  to  computer  frame 

4 

4>x 

X-component  of  vector  angle  from  true  to  platform  frame 

5 

(by 

Y-component  of  vector  angle  from  true  to  platform  frame 

6 

Z-component  of  vector  angle  from  true  to  platform  frame 

7 

5V, 

X-component  of  error  in  computed  velocity 

8 

Y-component  of  error  in  computed  velocity 

9 

Z-component  of  error  in  computed  velocity 

5h 

Error  in  vehicle  altitude  above  reference  ellipsoid 

11 

Shn 

Total  barometric  altimeter  correlated  error 

12 

User  clock  bias 

13 

^^Uclk 

User  clock  drift 
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Appendix  C.  Truth  Model  System  Dynamics  and  Process  Noise  Matrices 

The  LN-93  [26]  error-state  system  dynamics  matrix,  F,  is  a  93  x  93  matrix  with 
numerous  zero  elements.  This  matrix  is  parsed  out  into  portions  (submatrices)  that 
contain  non-zero  elements,  see  Equations  (64)  and  (65).  The  submatrices  of  the  reduced- 
order  truth  matrix  F  develop  by  Negast  [38]  are  shown  in  Tables  40  —  44.  A  complete 
listing  of  the  93-state  system  dynamics  matrix  can  be  found  in  the  Litton’ s  reference  [26] 
as  well  in  theses  at  AFTT  [38, 56].  The  system  dynamics  model  for  the  single-positioning 
GPS  and  DGPS  models  are  contained  within  Section  3.5. 

The  Litton  LN-93  truth  model  also  includes  a  93  x  93  process  noise  matrix,  Q  [26]. 
As  shown  in  Equations  (64)  and  (65),  the  WGN  process  noise  vector,  w,  was  subdivided 
into  subvectors,  Wi  and  W2,  corresponding  to  the  error  state  subvectors,  5xi  and  5x2. 
These  noise  subvectors  have  process  noise  covariance  (strength)  matrices  described  by 
Qii  and  Q22  listed  in  Tables  45  and  46.  The  process  noise  submatrices  for  the  GPS  truth 
models  are  contained  within  Section  3.5. 
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Table  40.  Elements  of  the  System  Dynamics  Submatrix,  Fn 


Element 

Term 

Element 

Term 

(1,3) 

-Py 

(1,8) 

~^RY 

IKSH 

Px 

mam 

^RX 

Py 

nfeH 

-P. 

-Q, 

(4,3) 

a 

y 

Oiin, 

(4,6) 

-0)in, 

(4,8) 

(5,1) 

Q, 

MEESM 

m. 

MHIM 

c 

^RX 

(6,1) 

-^y 

(6,2) 

(6,4) 

(6,5) 

-2Vp^-2V^a^ 

2F,Q. 

2Vfi, 

Him 

-4 

(7,7) 

‘~^z^RX 

(7,9) 

-Py-^^y 

(8,1) 

2VQ. 

X  y 

(8,2) 

-2V,a,-2Vfi, 

(8,3) 

2VJi, 

(8,6) 

-4 

Hn 

-2a, 

(8,8) 

~YzCrY 

p,+2a. 

(9,1) 

2VM, 

MIEM 

2V,G, 

IKBH 

(9,4) 

4 

(9,7) 

Py  +  2Qy  ^Yx^RX 

(9,8) 

(9,10) 

■dItM 

-K 

(9,12) 

-1 

*2 

(10,9) 

1 

mmm 

-K 

■dliMifli 

jtj-l 

1 

-1 

k. 

-hy 

(16,10) 

K 

(16,14) 

-K 

(16,16) 

K-\ 

=  Components  of  angular  rate,  navigation  frame  to  ECEF  frame 
^x,y,z  =  Components  of  angular  rate,  ECEF  to  inertial  frame 

ft) .  =  Components  of  angular  rate,  navigation  frame  to  inertial  frame 

^  =  Components  of  vehicle  velocity  vector  in  ECEF  coordinates 

^x,y,z  =  Components  of  specific  force  in  sensor  frame 

kjaM  -  Vertical  channel  gains 
a  =  Equatorial  radius  of  the  earth  (20,926,470  ft) 

Crx,ry  =  Components  of  earth  spheroid  inverse  radii  of  curvature 
^0  =  Equatorial  gravity  (32.08744  ft/sec^) 
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Table  41.  Elements  of  the  System  Dynamics  Submatrix,  F12 


Element 

Term 

Element 

Term 

Element 

Term 

(7,17) 

Cn 

(7,18) 

C\2 

Ci3 

(7,20) 

1 

(8,17) 

C2I 

(8,18) 

C22 

C23 

(8,21) 

1 

Csi 

mmm 

C32 

(9,19) 

C33 

1 

(9,11) 

k2 

(10,11) 

h 

(15,11) 

-h 

■iMIM 

kiieoo 

Note:  For  the  above  element  definitions,  to  =  0 

C  =  Coordinate  transformation  matrix  from  body  frame  to  navigation  frame, 


Table  42.  Elements  of  the  Dynamics  Submatrix,  F13 


1  Element 

Term 

Element 

Term 

Element 

Term 

Ql 

(4,24) 

Q2 

Qs 

(4,27) 

mmm 

Qi 

C22 

Q3 

(5,26) 

^21^  in^ 

mmm 

C31 

(6,24) 

Q2 

C33 

(6,27) 

^33^  in. 
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QCi.  Qd.  ca.  CO. 


Table  43.  Elements  of  the  Dynamics  Submatrix,  F14 


Element 

Term 

Element 

Term 

Element 

Term 

(7,29) 

(7,30) 

Q2 

(7,31) 

Cj3 

(7,32) 

(7,33) 

Cn^' 

(7,34) 

CnAf 

(7,35) 

c„K' 

(7,36) 

Q2 

(7,37) 

C.3 

Af 

1 

(7,38) 

(7,39) 

-CnA! 

(7,40) 

QA' 

(7,41) 

CuAf 

(8,29) 

Qi 

(8,30) 

C22 

(8,31) 

Q3 

(8,32) 

QiAf 

(8,33) 

Q2A 

(8,34) 

QsAf 

(8,35) 

Qi 

Af 

(8,36) 

C22 

1 

(8,37) 

c  * 

'^23 1 

Af 

(8,38) 

QA' 

(8,39) 

Q2^x 

(8,40) 

Qsa;  ■ 

(8,41) 

Q3A," 

(9,29) 

Qi 

(9,30) 

Q2 

(9,31) 

Q3 

(9,32) 

r"  4® 

(9,33) 

C32A; 

(9,34) 

C33Af 

(9,35) 

r 

^31 

(9,36) 

c 

^^32 

K 

(9,37) 

r 

'"33 

Af 

(9,38) 

Qia; 

(9,39) 

— r' 

^32^jc 

(9,40) 

C33A; 

(9,41) 

Axy  z  ~  Components  of  acceleration  in  the  body  frame 

A®  =  Specific  force  component  (includes  gravity) 


Table  44.  Elements  of  the  Dynamics  Submatrix,  F22 


Element 

Term 

Element 

Term 

Element 

Term 

(17,17) 

-K 

(18,18) 

(19,19) 

(20,20) 

~^Sg, 

(21,21) 

-K 

(22,22) 

~^Sg, 

(11,11) 

~Psh, 

.  =  Gyro  inverse  correlation  time  constants  (5  min) 

^  =  Accelerometer  inverse  correlation  time  constants  (5  min) 

•>€  -V 

^  =  Gravity  vector  error  inverse  correlation  time  constants  (V/20NM) 

^  =  Barometer  inverse  correlation  time  (10  min) 
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Table  45.  Elements  of  Process  Noise  Submatrix  Qn 


Table  46.  Elements  of  Process  Noise  Submatrix  Q22 


=  PSD  value  of  gyro  drift  rate  white  noise  (6.25  x  10"*** 

=  PSD  value  of  accelerometer  white  noise  (1.037  x  10"^  -^) 

see 
1 3 

=  Variances  of  gyro  drift  correlated  noise  (3.086  x  10' 

=  Variances  of  accelerometer  correlated  noise  (4.147  x  10’^ 

see 

-‘6  2 

=  Variances  of  gravity  vector  error  component  correlated  noise  (1.93  x  10’  deg  ) 
=  Variance  of  barometer  correlated  noise  (10,000  ft^) 
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Appendix  D.  Filter  Model  Process  and  Measurement  Noise  Matrices 


This  appendix  provides  the  values  of  the  system  dynamics  driving  noise  employed  in 
the  M^AE’s  Kalman  filters  used  in  the  software  simulations.  These  values  are  adjusted 
from  the  levels  specified  in  the  truth  model.  This  is  done  to  compensate  for  the 
inadequacies  introduce  by  using  a  small  subset  of  states  to  represent  the  behavior  of  a 
higher-order  truth  model.  Tables  47  -  49  below  shows  the  diagonal  (autocovariance) 
terms  of  the  process  and  measurement  noise  matrices,  Q  and  R,  respectively.  Off- 
diagonal  terms  (cross-covariance)  are  assumed  to  be  zero. 


Table  47.  Process  Noise  Strength  Values  (Q)  for  Filter  States  (0.4  nmi/hr  INS) 


State 

DGPS 

Units  1 

5d, 

56, 

^bswm 

56, 

0.0 

20.0 

HIHHH 

28.0 

85.0 

wimmim 

5V, 

5Vy 

sv. 

33,000.0 

mmi 

16.0 

Shs 

2.0 

mmm 

0.1 
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Table  48.  Process  Noise  Strength  Values  (Q)  for  Filter  States  (4.0  nmi/tir  INS) 


State 

DGPS 

Units  1 

56, 

■bseh 

56, 

fTffMI 

56, 

0.0 

14,500.0 

14,500.0 

520.0 

5V, 

500.0 

SV, 

500.0 

hbh 

5V, 

43,000.0 

6h 

18.0 

IHSHH 

5hB 

10.0 

HSHH 

Table  49.  Measurement  Noise  Strengths  Values  (R)  for  the  Truth  and  Filter  States 


Measurement 

Truth  Model  Value 

Filter  Model  Value 

Units  1 

Barometric  Altimeter 

2,500.0 

3,500.0 

mm 

DGPS 

9.0 

RpRl2.3A 

mm 

GPS  Pseudolite 

1.0 

RpRS 

i^^nn 

Radar  Altimeter 

0.09 

0.12 

126 


Bibliography 


1.  Athans,  M.  and  C.  B.  Chang.  Adaptive  Estimation  and  Parameter  Identification 
Using  a  Multiple  Model  Estimation  Algorithm.  Technical  Note  1976-28:  Lincoln 
Laboratory,  Massachusetts  Institute  of  Technology,  Lexington,  MA,  June  1985. 
ESD-TR-76-184. 

2.  Bagley,  Daniel,  T.  GPS/INS  Integration  for  Improved  Aircraft  Attitude  Estimates. 

MS  Thesis,  AFTT/GE/ENG/91D-04,  School  of  Engineering,  Air  Force  Institute  of 
Technology,  Wright-Patterson  AFB,  OH,  December  1991. 

3.  Baram,  Yoram.  Information,  Consistent  Estimation  and  Dynamic  System 
Identification.  PhD  Dissertation,  Massachusetts  Institute  of  Technology,  Cambridge, 
MA,  November  1976. 

4.  Bohenek,  Brian  J.  The  Enhanced  Performance  of  an  Integrated  Navigation  System  in 
a  Highly  Dynamic  Environment.  MS  Thesis,  AFIT/GE/ENG/94D-01,  School  of 
Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH, 
December  1994. 

5.  Braff,  Ronald,  Description  of  the  FAA’s  Local  Area  Augmentation  System  (LAAS). 
NAVIGATION,  Vol.  44,  No.  4,  Winter  1997-98,  pp.  411-424. 

6.  Britting,  Kenneth  R.  Inertial  Navigation  Systems  Analysis.  New  York:  Wiley- 
Interscience,  1971. 

7.  Britton,  Ryan  L.  A  Differential  GPS  Aided  INS  for  Aircraft  Landings.  MS  Thesis, 
AFrr/GE/ENG/95D-03,  School  of  Engineering,  Air  Force  Institute  of  Technology, 
Wright-Patterson  AFB,  OH,  December  1995. 

8.  Cheng,  David  K.  Field  and  Wave  Electromagnetics,  2’'^  Ed.  Reading,  MA: 
Addison-Wesley  Publishing,  Inc.,  1992. 

9.  Defense  Mapping  Agency.  Department  of  Defense  World  Geodetic  System  1984. 
DMA  Technical  Report,  DMA  TR  8350.2,  30  September  1987. 

10.  Department  of  Transportation.  Automatic  Landing  Systems,  Federal  Aviation 
Administration,  AC  20-57 A,  Washington,  DC,  12  January  1971. 

11.  Department  of  Transportation.  Engineering  and  Operational  Issues  Associated  with 
the  Application  of  Satellite-Based  Navigation  to  Precision  Approach  and  Landing, 
Revision  A,  NAS  System  Engineering  Service,  Federal  Aviation  Administration, 
Washington,  DC,  February  1992. 


127 


12.  Department  of  Transportation.  FAR-AIM  (Federal  Aviation  Regulations  and 
Airman’s  Information  Manual)  Part  1,  Aviation  Supplies  &  Academics,  Inc.,  Renton, 
WA,  1993. 

13.  Gray,  Robert  A.  An  Integrated  GPS/INS/BARO  and  Radar  Altimeter  System  for 
Aircraft  Precision  Approach  Landings.  MS  Thesis,  AFrr/GE/ENG/94D-13,  School 
of  Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH, 
December  1994. 

14.  Grejner-Brzezinska,  Dorota  A.  “Direct  Platform  Orientation  with  Tightly  Integrated 
GPS/INS  in  Airborne  Applications, "  1998. 

15.  Griffin,  Gordon  C.  and  Peter  S.  Maybeck.  “MMAE/MMAC  Techniques  Applied  to 
Large  Space  Structure  Bending  with  Multiple  Uncertain  Parameters,”  Proceedings  of 
the  34'^  IEEE  Conference  on  Decision  and  Control,  New  Orleans,  LA,  1 153-1158, 
December  1995. 

16.  Gustafson,  John  A.  and  Peter  S.  Maybeck.  “Flexible  Spacestructure  Control  Via 
Moving-Bank  Multiple  Model  Algorithms,”  IEEE  Transactions  on  Aerospace  and 
Electronic  Systems,  Vol.  30,  No.  3, 750-757,  July  1994. 

17.  Hansen,  Neil  P.  Incorporation  of  Carrier-Phase  Global  Positioning  System 
Measurements  into  the  Navigation  Reference  System  for  Improved  Performance.  MS 
Thesis,  AF1T/GE/ENG/93D-40,  School  of  Engineering,  Air  Force  Institute  of 
Technology,  Wright-Patterson  AFB,  OH,  December  1993. 

18.  Hartman,  Randolph  and  Don  Johnson.  Demonstration  of  a  P(Y)-Code  Differential 
GPS  Precision  Approach  System.  NAVIGATION,  Vol.  45,  No.  1,  Spring  1998,  pp. 
39-50. 


19.  Himing,  James  L.  Optimal  Kalman  Filter  Integration  of  a  Global  Positioning  System 
Receiver  and  an  LN-94  Inertial  Navigation  System.  MS  Thesis,  AF1T/GE/ENG/90S- 
02,  School  of  Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB, 
OH,  September,  1990. 

20.  Hong,  Lang.  Professor  of  Electrical  Engineering,  Wright  State  University,  Fairborn, 
OH,  Course  Notes  from  EE-71 7  Multisensor  Data  Fusion.  June  1998. 

21.  Huangqi,  Sun,  et  al.  “An  Investigation  of  Airborne  GPS/INS  for  High  Accuracy 
Position  and  Velocity  Determination,”  Proceedings  of  the  Institute  of  Navigation 
1994  National  Technical  Meeting,  San  Diego,  CA,  801-809,  24-26  January  1994. 

22.  Jacob,  Thomas.  “Landing  System  Using  GPS/IMU  System  Integration,”  Proceedings 
of  the  Institute  of  Navigation  1994  National  Technical  Meeting,  San  Diego,  CA,  298- 
305,  24-26  January  1994. 


128 


23.  Johnson,  Gregory  B.  Closed-Loop  Performance  of  GPS  Aided  INS.  MS  Thesis, 
AFIT/GE/ENG/89D-19,  School  of  Engineering,  Air  Force  Institute  of  Technology, 
Wright-Patterson  AFB,  OH,  December  1989. 

24.  Kalman,  R.  E.  “A  New  Approach  to  Linear  Filtering  and  Prediction  Problems,” 
Transactions  ASME,  Series  D:  Journal  of  Basic  Engineering,  S2:34-35, 1960. 

25.  Kaplan,  Elliott  D.  Understanding  GPS,  Principles  and  Applications.  Boston:  Artech 
House,  Inc.,  1996. 

26.  Knudsen,  L.  Performance  Accuracy  (Truth  Model/Error  Budget)  Analysis  for  the 
LN-93  Inertial  Navigation  System  Inertial  Navigation  Unit.  Technical  Report:  Litton 
Guidance  and  Control  Systems,  January  1985.  Doc.  No.  469414,  DID  No.  Di-S- 
21433  B/T:  CDRLNo.  1002. 

27.  Lewantowicz,  Zdzislaw  H.  and  Randall  N.  Paschall.  '‘Deep  Integration  of  GPS,  INS, 
SAR,  and  Other  Sensor  Information,”  1995. 

28.  Lund,  Elvind  J.  On-line  Discrimination  and  Estimation  in  Multiple  Regime  Systems. 
PhD  Dissertation,  Division  of  Engineering  Cybernetics,  The  Norwegian  Institute  of 
Technology,  University  of  Trondheim,  Trondheim,  Norway,  June  1992. 

29.  Lund,  Elvind  J.,  et  al.  “Multiple  Model  Estimation  with  Inter-Residual  Distance 
Feedback,”  Mode/mg,  Identification  and  Control,  Vol.  13,  No.  3, 127-140, 1992. 

30.  Maybeck,  Peter  S.  Stochastic  Models,  Estimation  and  Control,  I.  New  York: 
Academic  Press,  Inc.,  1982. 

31.  Maybeck,  Peter  S.  Stochastic  Models,  Estimation  and  Control,  II.  New  York: 
Academic  Press,  Inc.,  1982. 

32.  Maybeck,  Peter  S.  Stochastic  Models,  Estimation  and  Control,  III.  New  York: 
Academic  Press,  Inc.,  1982. 

33.  Miller,  Mikel  M.  Modified  Multiple  Model  Adaptive  Estimation  (M^AE)  for 
Simultaneous  Parameter  and  State  Estimation.  PhD  Dissertation,  AFrr/GE/ENG/98- 
02,  School  of  Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB, 
OH,  March  1998. 

34.  Mosle,  William  B.  Detection,  Isolation,  and  Recovery  of  Failures  in  an  Integrated 
Navigation  System.  MS  Thesis,  AFTr/GE/ENG/93D-28,  School  of  Engineering,  Air 
Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH,  December  1993. 

35.  Musick,  Stanton  H.  PROFGEN  -  A  Computer  Program  for  Generating  Flight 
Profiles.  Technical  Report,  Air  Force  Avionics  Laboratory,  Wright-Patterson  AFB, 
OH,  November  1976.  AFAL-TR-76-247,  DTIC  ADA034993. 


129 


36.  Musick,  Stanton  H.  and  Neil  Carlson.  User's  Manual  for  a  Multimode  Simulation  for 
Optimal  Filter  Evaluation  (MSOFE).  Air  Force  Avionics  Laboratory,  Wright- 
Patterson  AFB,  OH,  April  1990.  AFWAL-TR-88-1138,  AFWA17AARN-2. 

37.  NAVSTAR  GPS  Joint  Program  Office  (SMC/CZ).  NAVSTAR  GPS  User  Equipment 
Introduction.  Space  and  Missile  Systems  Center,  Los  Angeles  AFB,  CA.  February 
1991. 

38.  Negast,  William  J.  Incorporation  of  Differential  Global  Positioning  System 
Measurements  Using  an  Extended  Kalman  Filter  for  Improved  Reference  System 
Performance.  MS  Thesis,  AFIT/GE/ENG/91D-41,  School  of  Engineering,  Air  Force 
Institute  of  Technology,  Wright-Patterson  AFB,  OH,  December  1991. 

39.  Nielsen,  Robert  L.  Development  of  a  Performance  Evaluation  Tool  (MMSOFE)  for 
Detection  of  Failures  with  Multiple  Model  Adaptive  Estimation  (MMAE).  MS  Thesis, 
Wright-Patterson  AFB,  OH,  December  1993. 

40.  Nisner,  Paul  and  John  Owen.  “Practical  Measurements  of  Radio  Frequency 
Interference  to  GPS  Receivers  and  an  Assessment  of  Interference  Levels  by  Flight 
Trials  in  the  European  Regions,”  Proceedings  of  Institute  Of  Navigation  GPS-95, 
Parti,  1373-1382, 12-15  September  1995. 

41.  Parkinson,  Bradford  W.  and  Per  K.  Enge.  “  Differential  GPS.”  Global  Positioning 
System:  Theory  and  Applications.,  Vol.  H,  Ch.  1,  Washington  DC:  American  Institute 
of  Aeronautics  and  Astronautics,  Inc.,  1996. 

42.  Parkinson,  Bradford  W.,  “GPS  Error  Analysis.”  Global  Positioning  System:  Theory 
and  Applications.,  Vol.  I,  Ch.  11,  Washington  DC:  American  Institute  of  Aeronautics 
and  Astronautics,  Inc.,  1996. 

43.  Pinson,  J.  C.  “Inertial  Guidance  for  Cruise  Vehicles,”  Guidance  and  Control  of 
Aerospace  Vehicles  (C.  T.  Leondes,  ed,).  McGraw-Hill,  NY,  1963. 

44.  Ritland,  John  T.  “Impact  of  Inertial  System  Quality  on  GPS-Inertial  Performance  in  a 
Jamming  Environment,”  Proceedings  of  the  American  Institute  of  Aeronautics  and 
Astronautics,  1459-1467, 1987. 

45.  Rowson,  Stephen  V.,  et  al.  “Performance  of  Category  IIIB  Automatic  Landings 
Using  C/A  Code  Tracking  Differential  GPS,”  Proceedings  of  the  Institute  of 
Navigation  1994  National  Technical  Meeting,  San  Diego,  CA,  759-767, 24-26 
January  1994. 

46.  Schiller,  Gregory  J.  and  Peter  S.  Maybeck.  “Space  Structure  Control  Using  Multiple 
Model  Adaptive  Estimation  and  Control,”  Proceedings  of  the  13^^  IF  AC  Symposium 
Automatic  Control  in  Aerospace-Aerospace  Control  ‘94,  September  1994. 


130 


47.  Schnedorf,  Jerome  G.,  Ill  Components  of  the  Instrument  Landing  System.  World 
Wide  Web,  http://home.sprvnet.com/sprvnet/iavschne/Ils2.htm.  1997. 

48.  Sheldon,  Stuart  N.  An  Optimal  Parameter  Discretization  Strategy  for  Multiple  Model 
Adaptive  Estimation  and  Control.  PhD  Dissertation,  School  of  Engineering,  Air 
Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH,  December  1989. 

49.  Sheldon,  Stuart  N.  and  Peter  S.  Maybeck.  “An  Optimizing  Design  Strategy  for 
Multiple  Model  Adaptive  Estimation  and  Control,”  IEEE  Transactions  on  Automatic 
Control,  Vol.  38,  No.  4, 651-654,  April  1993. 

50.  Sokol,  Charles  W.  Performance  Evaluation  of  a  GPS-Aided  INS  in  a  Filter-Driving- 
Filter  Configuration.  MS  Thesis,  AFrr/GE/ENG/94M-02,  Air  Force  Institute  of 
Technology,  Wright-Patterson  AFB,  OH,  March  1994. 

51.  Stacey,  Richard  D.  A  Navigation  Reference  System  (NRS)  Using  Global  Positioning 
System  (GPS)  and  Transponder  Aiding.  MS  Thesis,  AFlT/GE/ENG/9  lM-04,  School 
of  Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH, 

March  1991. 

52.  The  MathWorks,  Inc.,  21  Elliot  Street,  Natick,  MA  01760.  MATLAB.  December 
1997.  Version  5.2. 

53.  Vanek,  Barry  J.  GPS  Signal  Offset  Detection  and  Noise  Strength  Estimation  in  a 
Parallel  Kalman  Filter  Algorithm.  MS  Thesis,  AFIT/GE/ENG/99M-30,  School  of 
Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH,  In 
Progress,  March  1999. 

54.  Vasquez,  Juan  R.  Detection  of  Spoofing,  Jamming,  or  Failure  of  a  Global 
Positioning  System  (GPS).  MS  Thesis,  AFrr/GE/ENG/92D-37,  School  of 
Engineering,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH, 
December  1992. 

55.  Vasquez,  Juan  R.  New  Algorithms  for  Moving-Bank  Multiple  Model  Adaptive 
Estimation.  PhD  Dissertation,  AFIT/DS/ENG/98-10,  School  of  Engineering,  Air 
Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH,  May  1998. 

56.  White,  Nathan  A.  MMAE  Detection  of  Interference/ Jamming  and  Spoofing  in  a 
DGPS-Aided  INS.  MS  Thesis,  AFrr/GE/ENG/96D-21,  School  of  Engineering,  Air 
Force  Institute  of  Technology,  Wright-Patterson  AFB,  OH,  December  1996. 

57.  Young,  Brian  J.  An  Integrated  Synthetic  Aperture  Radar /Global  Positioning  System 
/Inertial  Navigation  System  for  Target  Geolocation  Improvement.  MS  Thesis, 
AFrr/GE/ENG/99M-32,  School  of  Engineering,  Air  Force  Institute  of  Technology, 
Wright-Patterson  AFB,  OH,  In  Progress,  March  1999. 


131 


58.  Zumwalt,  Michel  P.,  Maj.  Test  Pilot  Candidate,  Edwards  AFB,  CA.  Personal 
Interview.  October  1998. 


132 


Vita 


Steve  Chastain’s  hometown  is  Orlando  FL.  He  received  a  Bachelor  of  Science  degree 
in  Electrical  Engineering  from  Florida  Atlantic  University.  He  entered  Officer  Training 
School  at  Lackland  AFB  TX  and  received  his  commission  on  22  September  1993.  His 
first  Air  Force  assignment  was  with  the  Space-Based  Infrared  System  Program  Office  at 
Los  Angeles  AFB  CA.  His  duties  included  overseeing  the  testing  and  maintenance  of  the 
infrared  sensors  on-board  the  Defense  Support  Program  early-warning  satellites.  He  then 
spent  a  year  matrixed  to  the  Aerospace  Corporation  in  the  Astrodynamics  Department 
doing  studies  on  the  atmospheric  effects  in  low  earth  orbit  and  designing  a  orbital 
simulator  for  molniya  (highly  eccentric)  orbits.  In  August  1997,  he  assigned  to  attend  the 
Guidance  &  Control  program  at  the  Air  Force  Institute  of  Technology  at  Wright- 
Patterson  AFB  OH,  specializing  in  stochastic  modeling  and  navigation  systems.  After 
graduating  on  23  March  1999,  his  follow-on  assignment  was  to  the  746**'  Test  Squadron 
at  Holloman  AFB  NM. 


133 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  Information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 


3.  REPORT  TYPE  AND  DATES  COVERED 


1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE 

March  1999 


4.  TITLE  AND  SUBTITLE 

Integrated  GPS/INS  Precision  Approach  Landing  with  M3AE  Interference  Avoidance 


6.  AUTHORCS) 
Stephan  H.  Chastain 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 
Air  Force  Institute  of  Technology 
2950  P  Street 

Wright-Patterson  AFB  OH  45433-7765 


MS  Thesis 


5.  FUNDING  NUMBERS 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFIT/GE/ENG/99M-03 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 
Juan  Vasquez,  Cpt,  USAF 
AFRL/SNAR  Bldg  620 
2241  Avionics  Circle 
Wright-Patterson  AFB  OH  45433-7321 

937-255-5668  x4014 _ _ 


1 1 .  SUPPLEMENTARY  NOTES 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 
Distribution  Unlimited 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

Several  past  research  efforts  have  developed  methods  to  take  advantage  of  Global  Positioning  System  (GPS)  positioning  and 
apply  it  to  a  precision  landing  system  (PLS).  There  have  been  proposals  to  phase  out  the  current  Instrument  Landing  System 
(ILS)  in  favor  of  a  more  cost-efficient  and  effective  system.  Accomplishments  have  been  made  in  detailing  a  system 
implementing  an  INS  aided  with  differential  GPS,  a  GPS  pseudolite,  and  a  radar  altimeter  to  handle  the  critical  PLS 
requirements.  This  research  applies  the  newly  developed  Modified  Multiple  Model  Adaptive  Estimation  (M3AE)  architecture 
in  an  attempt  to  enhance  the  performance  of  a  PLS  in  an  environment  involving  GPS  interference.  The  M3AE  uses  Multiple 
Model  Adaptive  Estimation  (MMAE)  and  Kalman  filtering  techniques  to  estimate  the  levels  of  interference  and  the  navigation 
performance  of  the  aircraft  simultaneously.  In  addition,  in  the  original  development  of  M3AE,  the  truth  and  filter  model  used 
were  of  the  same  order.  This  research  serves  as  a  demonstration  of  M3AE  applied  to  system  where  the  truth  model  is  of  a 
higher  order  than  the  filter  model. 


14.  SUBJECT  TERMS 


Precision  Approach  Landing;  Inertial  Navigation  Systems  (INS);  Global  Positioning  System 
(GPS);  Kalman  Filtering;  Adaptive  Filtering 


15.  NUMBER  OF  PAGES 
146 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRAC 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


Unclassified 


Unclassified 


Unclassified 


Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


